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Abstract: The axino, the fermionic superpartner of the axion, is a well-motivated can- 
didate for cold dark matter if it is the lightest supersymmetric particle. Since the axino 
couples very weakly to the matter multiplets, the next-to-lightest supersymmetric particle 
(NLSP) has a long lifetime, which has important consequences for both cosmology and 
collider phenomenology. Assuming that a charged slepton is the NLSP, we calculate the 
complete leading one- and two-loop contributions to its decay. We analyze in detail con- 
straints on the parameters space from cosmology and discuss how this scenario can be 
probed at colliders. Scenarios in which both the axino and the gravitino are lighter than 
the long-lived charged slepton are also explored with particular emphasis on cosmological 
constraints and collider phenomenology. 
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1. Introduction 

There are compelling reasons for physics beyond the standard model (SM). As the exper- 
iments at the Large Hadron Collider (LHC) arc taking data, we could well see the first 
hints of this new physics soon. Along with high energy experiments, cosmology provides 
important testing grounds for models of new physics and their effects in the early universe. 
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Among the extensions of the SM, supersymmetry (SUSY) is in some sense the most 
ambitious and appeahng one [1-3]. It evokes a completely new symmetry of nature and 
opens new angles on old puzzles, such as the hierarchy problem, the unification of the 
fundamental forces or even the incorporation of gravity into the SM of particle physics. 

Supersymmetry is particularly attractive in connection with the Peccei-Quinn (PQ) 
mechanism [4,5]. Relying on a global U(1)pq symmetry that is broken spontaneously 
at a high scale /a, this mechanism provides a compelling potential solution of the strong 
CP problem, which consists of the question of why CP violation in QCD is so small. 
The breaking of the PQ symmetry gives rise to the famous axion which could provide 
(part of) dark matter [6-8] or influence other physics because of its lightness and its weak 
interactions. In fact, axion searches and a large variety of astrophysical and cosmological 
considerations point to the following lower limit on the PQ scale [9, 10]: 

/„>6xlO«GeV. (1.1) 

Because of the large breaking scale, new physics must be involved in viable axion 
models and the coupling of the axion to the SM fields must be highly suppressed. One 
class of such models are 'hadronic models' - also known as KSVZ axion models [11,12] - 
where new heavy 'quark' fields are introduced that couple to the axion. Via interactions 
of gluons (and possibly also of other gauge fields) with those heavy coloured fields, loop- 
suppressed axion couplings to the SM fields emerge. SUSY versions of this type of models 
are the ones on which we focus in this work. 

In the considered SUSY theories, there must be a fermionic partner to the axion: the 
axino [13-15].^ Its mass depends on the details of the model [13,20-27] and can be such 
that the axino is the lightest super symmetric particle (LSP), which is assumed in most of 
this work together with R-parity conservation. In such settings, the axino can constitute 
a (significant) fraction of dark matter [22,28-36]. 

With the axino being lighter than the SM superpartners, the lightest sparticle in the 
minimal supersymmetric standard model (MSSM) is unstable and can thus be a charged 
slepton such as the lighter stau [32,35,37-39]. Such scenarios are particularly attractive 
since the charged slepton would typically be long-lived. Its existence could thus be asso- 
ciated with striking signatures of tracks of heavy quasi-stable charged particles in collider 
experiments. Some fraction of the charged sleptons are even expected to get stopped 
in the main detectors or the surrounding material, such that their decays can be anal- 
ysed [37, 40-44]. Long-lived charged particles can also affect big bang nucleosynthesis 
(BBN) since the decay products can reprocess primordial nuclei [39]. Moreover, being elec- 
trically charged, they themselves can form primordial bound states, which may catalyse 
the synthesis of ^Li and ^Be in the early universe [35,45-47]. 

Given the theoretical as well as the experimental attractiveness of axino dark matter 
scenarios with a long-lived charged slepton, we present in this paper a comprehensive study 

^ There is also the scalar partner of the (pseudoscalar) axion, the saxion, whose late decays could lead 

to significant non-standard entropy production [16-19]. For the present investigation, however, we assume 
a standard thermal history and thereby implicitly a sax;ion mass such that sax;ion eflFects are negligible. 
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of this class of models. In Sec. 2 we describe the considered SUSY hadronic axion models. 
They typically have two scales, the weak scale and the breaking scale of the PQ symmetry, 
which differ by several orders of magnitude. This leads to possibly large logarithms at 
the two-loop level which require detailed analyses. In Sec. 3 we therefore present a two- 
loop calculation of the decay rate of the long-lived charged slepton and an analysis of the 
various relevant decays. This is necessary to correctly interpret cosmological constraints 
and collider studies of axino signatures and the PQ scale /a. While some of our results 
are presented in previous short communications [35,39], we provide here a more complete 
systematic treatment along with details of the calculations. We also note that the applied 
methods will be useful for other two-scale models and we have thus assembled the technical 
details in the appendix. 

In Sec. 4 we outline the considered cosmological settings and analyse the associated 
cosmological constraints. These include constraints from the present dark matter density, 
from structure formation and from primordial nucleosynthesis. The phenomenology at 
the LHC and at a future linear collider such as the International Linear Collider (ILC) is 
explored in Sec. 5. Here we propose and study four benchmark scenarios in the constrained 
MSSM (CMSSM) extended with a hadronic axion model. Moreover, we address prospects 
for analysing late charged slepton decays, for the microscopic measurements of the PQ 
scale /a and for distinguishing between the axino LSP and the gravitino LSP, which is 
another well- motivated dark matter candidate that allows for the existence of a long-lived 
charged slepton; cf. [48] and references therein. 

In Sec. 6 we study settings in which both the axino and the gravitino are lighter than 
the long-lived charged slepton, which is still considered to be the lightest SM superpartner. 
Such settings can occur naturally in gauge-mediated SUSY breaking scenarios [26,27]. Here 
we also provide results on partial decay widths and the lifetime of the long-lived charged 
slepton, on the associated cosmological constraints and on the possible impact on collider 
phenomenology. Conclusions are presented in the final section. 

2. Supersymmetric Hadronic Axion Models 

In line with our previous studies [35,39], we consider the KSVZ-type or hadronic axion 
models [11,12] in a SUSY setting; see also [14,15,26,27], or [13,49,50] for other supersym- 
metrised axion models. Here new heavy quark multiplets Qi^2 are introduced that couple 
to the axion multiplet Only the new (s)quarks Qi,2 and ^> arc charged under U(1)pq. 
The corresponding interactions are described by the following superpotential 

WpQ = y^QiQ2 (2.1) 

with the Yukawa coupling y. Interactions of MSSM fields - which do not carry PQ charge 
in KSVZ-type models - with the axion multiplet $ are mediated via the Qi^2 fields, which 
are coloured and may also carry other charges under the SM gauge group. 

While our study can easily be generalized, we focus on the case in which Qi,2 are 
SU(3)c triplets and SU(2)l singlets, with electrical charge eg and hypercharge ey = eg, as 
given in Table 1. With the 2-component field notation in Table 1, the 4-component fields 
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chiral multiplet 


U(1)PQ 


(SU(3)e,SU(2)L)Y 


$ = + + ^K^^ 


+ 1 


(1,1)0 


Qi = Qi + ^f2q^d + Fi^^ 


-1/2 


(3,l)+eg 




-1/2 


(3*,l)-eQ 



Table 1: The axion multiplet the heavy KSVZ quark multiplets (3i,2, and the associated 
quantum numbers considered in this work. The complex scalar </> = s + ia is a combination of the 
saxion s and the axion a. 

describing the axino and the heavy KSVZ quark arc given respectively by 

and Q = ( . (2.2) 

For the heavy KSVZ (s)quark masses, we use the SUSY limit 

Mq = Mq^^^ =y{^)=y fa/V2, (2.3) 

with both y and fa taken to be real by field redefinitions. The lower limit on the PQ 
scale (1.1) then implies, for y = 0(1), a large hierarchy between the mass of the new 
KSVZ fields and the weak and the soft SUSY mass scales: 

-^Q/Q, o ^ "^SUSY- (2.4) 



1,2 



For clarification of the definition of fa = \/2(0) in the considered models, note that 
axion-gluon, axion-photon, axino-gluino-gluon and axino-bino-photon/Z boson interac- 
tions are obtained as described by the effective Lagrangians 

^agg{g) = 32^2 ^ G^^C'^'' , (2.5) 



2 



^aBy/z = i [^M,^-] B {cosOwF^, - sinOwZ^,) , (2.8) 

once the heavy KSVZ (s) quarks are integrated out. Here G^,y, F^i, and Z^j^ are the field 
strength tensors of the gluon, photon and Z boson, respectively, with the duals = 
^ij.upaG"''"^ and -F^j, = e^^p^PP'^ /2, the strong coupling gs, e = y/Aira, and the Weinberg 
angle 9w- The weak hypercharge of the heavy KSVZ fields eq determines CaVY = Gcg 
and Ca-y-y, which becomes Cayy = Gcq — 2(4 + z)/[3(l + z)] after chiral symmetry breaking, 
where z = rriu/md ~ 0.56 denotes the ratio of the up and down quark masses. 

While the axino mass depends on the details of the model [13,20-27], we simply 
treat it as a free parameter. When focussing on the axino LSP case, this parameter is 
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bounded from above by the mass of lightest MSSM superpartner, which is considered to be 
the hghtest charged slepton in all of this work. For simplicity and concreteness, we focus 
mainly on scenarios where this slepton is a dominantly right-handed lighter stau: ri « tr.^ 
This is often a good approximation at least for small values of the mixing angle tan /3 in 
the Higgs sector (which is the ratio of the two MSSM Higgs doublet vacuum expectation 
values). Moreover, while the weak scale mass hierarchy = m^^ < ~ m^j^ is typically 
predicted in scenarios such as the CMSSM, most of our results can be directly transferred to 
the alternative cases in which the lightest slepton is the right-handed selectron cr or smuon 
juR. It will also be straightforward but require more effort to generalize the calculations 
of the decays to cases in which the lightest slepton Zi has a left-handed component Zl, 
including the case of a lighter stau t\ = cos tl -|- sin ^f^Tl with cos Of > 0. 

For the h = Zr cases, the coupling of the respective lighter charged slepton to neu- 
tralinos is dominated by the bino coupling. Here we assume for further simplicity that the 
mixing in the neutralino sector is such that one of the neutralino states is an (almost) pure 
bino. In fact, the cosmological considerations and collider studies are presented for spectra 
in which that state is the lightest neutralino, Xi = while the results presented in Sees. 3 
and 6.1 are not restricted to this case. 

If SUSY is realised as a gauge symmetry [51], also the gravitino G - the superpartner 
of the graviton - will exist as the gauge field associated with local SUSY transformations. 
Its mass itIq depends on the SUSY breaking mechanism and can thus often be treated as 
a free parameter. For example, rriQ can be in the eV range in gauge-mediated breaking 
schemes [52-54] and in the GeV to TeV range in gravity-mediated breaking schemes [55-57] . 
Accordingly, the following mass orderings are possible: (i) rria < mj^ < uIq, (ii) uIq < 
rrii^ < nia, (iii) rria < < m^-^ and (iv) < nia < mj^. In most of this work we will 
assume ordering (i). Here additional cosmological constraints may occur since gravitinos 
produced thermally in the early universe [58] may decay during/after BBN [59-61]. In 
Sec. 5.2, case (ii) is considered when exploring collider prospects to distinguish between 
the axino LSP and the gravitino LSP. Moreover, a dedicated treatment of the cases (iii) 
and (iv) is presented in Sec. 6, where we will also study cosmological constraints associated 
with the gravitino. For those gravitino-related constraints in the cases (i) and (ii), we refer, 
e.g., to Refs. [59-61] and [47], respectively. 

3. Decays of Chcirged NLSP Sleptons into Axino Dcirk Matter 

In this section, we outline the calculations and present the results of the 2-, 3- and 4-body 
decays of the stau NLSP into the axino. In the particle physics setting considered in this 

paper, this is the 2-body decay of the right-handed stau into the tau and the axino, tr ra, 
the 3-body decay into the tau, the axino and a photon, tr — t- tu'J, and the 4-body decay 
into the tau, the axino and a qq pair, tr — t- raqq. Finally, we look at the resulting lifetime 
of the stau NLSP, which is governed by the 2-body decay. After the obvious substitutions, 

^When considering prospects for collider phenomenology in Sec. 5.1, CMSSM benchmark scenarios are 
considered in which the lighter stau n has a left-handed component of < 10%. 
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the results presented in this section axe equahy vahd for the alternative cr and /xr NLSP 
cases. 



3.1 2-Body Decays of Charged NLSP Sleptons 

Let us start with the 2-body decay tr ra. The relevant Feynman diagrams are shown 
in Fig. L As <C m^, we work in the limit rrir — >■ 0. 

For the calculation of the two-loop diagrams, we can make use of the large hierar- 
chy (2.4) between the mass of the heavy quark multiplot and the weak scale or the SUSY 
breaking scale, see (2.3). As a result, one can perform a heavy mass expansion (HME) 
in powers of 1/fa- Already the leading term oc 1/ fa of the amplitude provides a precise 
approximation since the sub-leading terms are suppressed by many orders of magnitude. 
After the HME, only vacuum two-loop integrals (i.e. integrals with zero external momenta) 
and one-loop integrals remain in the expansion coefficients, for which well-known results 
are available in literature. More details about the HME of the two-loop diagrams can be 
found in Appendix A. The methodology was laid out previously in [62,63]. 

After a lengthy calculation, we obtain the following result for the decay rate 



r(TR -> ra) 



rriril - Aa) 



IGtt 



(3.1) 



A 



STT^ctv fa 



3^ Ab loe 



2 ml 



'A-n 



, 2(1 - A^)[^a{l - A^) - ^(1 - Ag)] A^ - 1 



l-A-a 



+ 



+ 



2cl,^{Aa - A^){1 - y^A~/Aa) - A-a 

la ^B 



1 - Aa "° A^ 

2y/Aa{l - A^Y - ^{2 - A^){1 - A-a) + 2c%j{l - A^){Aa^ + A^y/A^a) 



{l-Aa){l-A^) 



log Ai 



Sw(V^ \/^b)[2(V^+ ^/A^Y Az{1 + y/AaA^)^ ^2 ^ _2 n ™2 _2 _2 



2 ri / Z Z n 2 2 Z \ 



I -An. 



m~ Co{m~, m~, 0, m~, 0, m~) 



2c^{VA^+^){Aa-A^) 
l-Aa 

2s^[-\/3a(l — Ag) - y^A^{l — Aa)] •! ^ (n. '1 2 n 2 2 



\~Aa 



■ ruz Co(0, m~, m~, 0, m^, m~) 



+ 



Sy\/VA^i2-^ AaA^ + 2A^ Az) („ , 2 2 2\ 13 /n n 2^^ 

iZrj: [Bo{m^, rnz, m~) - Bo{0, 0, m~) j 

*w(2\/^B "I" ^VAa — VAqAz) / „.2 2 2\ r, n 2\\ 

^ [Bo{m~, mi,m$) - Bo{0, 0, m$)j 

2[sl,VA'a{l - A^) - (4 + cl,)^{l - Aa)] 



1-An 



(3.2) 
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Figure 1: Feynman diagrams for the 2-body stau NLSP decay tr — >■ ra in the considered type of 
KSVZ axion models. The thick sohd (dashed) lines indicate heavy (s)quark propagators. 



with 



Aa — — f ) A^ — Az — ^, (3-3) 



m~ m~ m~ 

T T T 



cw = cos^w and sw = sin^w Here Bq and Co are the usual scalar one-loop self-energy 
and vertex functions, respectively [64]. We have arranged the result so as to make the 
independence of the regularisation scale (which appears implicitly in the Bq functions) 
manifest. 

A striking feature of (3.1) is the appearance of the logarithm term oc log[y^/^/(2m~)] 
in the first line of (3.2). For typical values, y ~ 1, /a ~ 10^ . . . 10^^ GeV, and ^ 
100 . . . 1000 GeV, the logarithm becomes large, log(y^/^/2m|) ~ 26 ... 63, and dominates 
the numerical result. We will further elaborate on this aspect below, and in particular in 
Sec. 3.4. 

Before proceeding, let us mention the sensitivity of r(TR — > r a) to the two large scales 
fa and ^QiQ^ 2' Depending on y, those scales may differ; cf. (2.3). This is an interesting 
property of the TRra vertex governed by the two-loop diagrams shown in Fig. 1. From 
the one-loop diagrams leading to the effective vertices described by (2.5)-(2.8), only a 
sensitivity to the scale fa appears. Phenomenological studies of both types of vertices may 
thereby probe the coupling y that relates the scales /„ and Mqiq^ ^ as can be seen in (2.3). 

3.2 3-Body Decays of Charged NLSP Sleptons 

Now we turn to the 3-body decay tr — )■ ra^y. As before, the tau mass will be neglected 
throughout. The 3-body decay occurs already at one-loop order, see Fig. 2 (a), where the 
heavy quarks and squarks run in the loop indicated by the blob. However, the two-loop 
contributions shown in Fig. 2 (b-j) can be numerically important since they could generate 
a large logarithmic factor log[y^/^/(2m~)], similar to the one in (3.1). This factor can 
partially compensate the loop suppression factor a. 

Again we calculate the associated matrix element of the diagrams using an asymptotic 
expansion for large (s)quark masses, (2.3), and retaining only the leading term oc 1/ fa- The 
complete result for the two-loop contribution is rather lengthy. However, since the 3-body 
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decay rate is suppressed with respect to that of the 2-body decay, it is sufficiently accurate 
to only keep the dominant logarithm term oc log[y^/^/(2m|)] of the two-loop part. The 
resulting differential decay rate including the one- and two-loop amplitudes reads 

d^r(?R ^ ra^) ^ _9_ e% x^jl - Ag - x^) 

dx^dcos^ 647r4 /2 [2 - - cos^)]2 ^ ' ^ 

"2x^(1 - Aa- Xj)[l -I- cos 6* Aa{l - COS e)] [1 + cose + A^{1 - cos e)] 
[ [2Aa - A^{2 - + x^cosO) + x^{l + cos6')]2 

^ 3a (llfl\ - Ag) + ^1^(1 -Ag- x^)]{l + cos 9) + A-gAgXyjl - cos 9) 

TTC^ °^\2m~/ 2Aa — A^(2 — x^ + x^cosO) + Xy{l + cos9) 

2(y'fa\ 2(1 - Ag)il -Ag- x^)il + cos^) + .Tgj2 - (1 - - cos^)] ' 

87r2c4 °^ V2"^|y ^ x2(l-cos^)(l-^a-x-,) 

where 

x-y = 2E^/m^ , (3.5) 
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is the photon energy Ey divided by its maximum value and 9 is the angle between the 
photon and tau in the rest frame of the stau. 

It is interesting to note that only the contributions from the diagrams in Fig. 2 (f j) 
are logarithmically enhanced. In these diagrams the photon is radiated either from the 
external tau or stau line or from a four-vertex, so that the loop integral is essentially the 
same as for the 2-body decay. If the photon is radiated off the heavy (s)quark loop, as in 
Fig. 2 (d,e), an additional heavy (s)quark propagator is introduced in this loop. In fact, the 
leading term in the HME of these diagrams is not 0{l/fa) but 0(1//^) and thus strongly 
suppressed. Finally, the absence of a large logarithm log(y^/^/2m~) in the diagrams 
Fig. 2 (b,c) can be understood as follows: a logarithmic dependence on the heavy (s)quark 
mass corresponds to an ultraviolet (UV) divergence if the (s)quark mass becomes infinite, 
^Q/Qi2 ~ yfa/V^ c»- When the (s)quark loop is integrated out. Fig. 2 (b,c) turn 
into one-loop diagrams with four propagators in the loop. One can see from simple power 
counting that these one-loop terms are UV finite so that the full two-loop contribution 
from Fig. 2 (b,c) cannot generate any large logarithms. 

In [37], the 2- and 3-body decay rates were calculated using an effective theory where 
heavy (s)quark loops were integrated out; cf. (2.8). As mentioned above, this method 
leads to a logarithmic divergence in the result. As in [32], that logarithmic divergence was 
regulated with the cutoff /„, and only the dominant contributions were kept. Moreover, 
the authors of [37] introduced a factor ^ and a mass scale m to parametrise the uncertainty 
associated with this cutoff procedure. Our direct two-loop result now allows us to make a 
direct connection with the parameters of the underlying model, yielding agreement between 
(3.4) and the 3-body result of [37] for ^ = 1, m = \/2m^ly and, as mentioned below (2.8), 
CaYY = Scq. The comparison of the 2-body results can be found in [35] and will also be 
addressed in Sec. 3.4. 

The term in the last line of (3.4) with the squared logarithm has an infrared (IR) 
divergence for .x^ and a collinear divergence for cos0 — > 1. These divergencies would 
be cancelled by the virtual three-loop correction to the 2-body decay channel. However, 
the calculation of this three-loop contribution is a formidable task and beyond the scope of 
this paper. Furthermore, as will be shown below, the branching ratio of the 3-body mode 
is small and has a minor effect on the total lifetime. Therefore, the three-loop contribution 
to the 2-body decay mode will have a similarly small effect and can be safely neglected. 
Nevertheless, our calculation of the 3-body decay mode provides a meaningful description 
of the energy and angular distributions of the photon (as long as and 1 — cos Q do not 
become too small), which can be used for analysing this decay at colliders, see Sees. 5 
and 6.3 below. 

Because of the soft and collinear divergencies, we consider the integrated 3-body rate 
with cuts on the scaled photon energy, x^ > a;™*, and on the photon angle, cos 6 < 1 — x^"*, 

r(5n^™7;<;xr).y^^^_ ico.e .l^,^/' . (3.6) 
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0.2 0.4 , 0.6 0.8 1 100 120 140 160 180 200 

■i:^"' mf[GeV] 

(a) (b) 

Figure 3: Branching ratio of the 3-body decay tr — ra7 as a function of (a) the soft and colhnear 
cut parameters x^^"* and a;™* and (b) the masses to? and rUg. In both plots the other parameters are 
set to nia = 10 GeV, fa ^ 10" GeV, leg] = 1/3, y = 1, as well as m? = 100 GeV, = 110 GeV 
for the left figure, and a;™* = x™* = 0.1 for the right figure. 



Figure 3 shows the branching ratio of the integrated 3-body decay rate, 



BR(rR 



Taj; x^j,"*; Xg*^*) 



(3.7) 



with 



r[^t = r(?R ^ ra) + r(?R ^ tSt; x™*; x^"*). (3.8) 

As can be seen from Fig. 3 (a), within reasonable ranges of x^^^ and x™*, the 3-body decay 
contributes only C(l%) to the total decay rate. The latter is therefore well approximated 
by the 2-body contribution. For = rn^-fm- — )■ 1, the branching ratio (3.7) grows rapidly 
as shown in Fig. 3 (b) for x'^^ = Xg"' = 0.1. This is because the dominant, logarithmically- 
enhanced part of r(rR — >■ ra) oc m^m?~ while r(TR — t- ra'y) oc m~. However, the branching 
ratio always stays below about 3.5% for nir < 5TeV and xj^,"* = Xg"* = 0.1. As long as 
m~/m~ <C 1, the influence of on the branching ratio is small. The product yfa appears 
only in the logarithm but determines the relative contribution of the one- and two-loop 
amplitudes to the 3-body decay rate. For the parameters in Fig. 3, both the one- and two- 
loop amplitudes are of roughly equal importance. Varying yfa by one order of magnitude 
up or down changes the 3-body branching ratio by 2-10%. Nevertheless, this factor is 
important when considering the prospects for distinguishing between the axino LSP and 
the gravitino LSP, as will be discussed in the last two paragraphs of Sec. 5.2 below. 
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3.3 4-Body Decays of Charged NLSP Sleptons 

The 4-body decay modes tr — t- ra'y^^ tr — t- Tal~^l~ and tr — t- 
rTiqq have a neghgible impact on the total stau decay rate 
since they are suppressed by an additional factor a compared 
to the 3-body decay. However, for > 100 s, hadronic energy 
injection [60, 61, 65-70] from tr raqq decays can affect 
the abundances of primordial light nuclei significantly, see 
Sees. 4.3 and 6.2. The evaluation of BBN constraints from 
this effect requires the calculation of the hadronic 4-body 
decays. 

The Feynman diagrams for tr — > raqq can be deduced 
from the diagrams in Fig. 2 by splitting the final state photon 
into a qq pair. In addition there are diagrams where this 
photon is replaced by a Z boson. 

The one-loop contribution, see Fig. 4, has been calcu- 
lated with the same HME techniques explained in the previous sections. For the two-loop 
amplitude, we again only retain the terms oc log(y^/^/2m~). Following the arguments 
of the previous section, no such logarithmic term is generated from diagrams where the 
(7 ~^ qq) 01 {Z qq) is emitted from inside the loops. One has to consider only the two- 
loop diagrams where the 7, Z is emitted from the external stau or tau lines or is attached 
via a quartic coupling to the external and internal stau and the internal 7, Z (see also Fig. 1 
of [39]). As a result, one can directly use the result for the 3-body decay without having 
to recompute the complete amplitude. As the only necessary modifications, the stau-7/Z 
or tau-j/Z splitting functions have to be added and the stau or tau momentum in the 
loop integrals moved off-shell. In our calculation we include the width of the Z boson by 
introducing a complex mass in the propagator, — t- — imzTz- 

The differential decay rate as a function of the qq invariant mass mqq is given by 

dr(TR — >■ raqq) 1 1 



Figure 4: One of the Feyn- 
man diagrams of the 4-body 
stau NLSP decay tr — >• 
raqq in KSVZ axion models. 
The blob indicates the heavy 
(s)quark loop. 



drugq 32m| (27r) 



dm-rv / dila / dQ\ 



(3.9) 



x|fca| Wvl \kg\ \MiTR Taqq)\'^, 

where m-j-v is the invariant mass of the final state r and the intermediate F = 7, Z boson, 
fta and ka are the solid angle and 3-momentum of the axino in the tr rest frame, fly and 
k'y are the solid angle and 3-momentum of the V = j,Z boson in the rV rest frame, and 
flq and kg are the solid angle and 3-momentum of the final state quark in the V = j,Z 
rest frame. 

One of the solid angle integrations, say f dQa. is trivial and simply corresponds to a 
redefinition of the coordinate system. One additional azimuthal angle integration is also 
trivial since the incoming stau carries no polarisation. The integration over the angle of 
the final state qq pair yields, for massless quarks, 

167r 



dn'^ [uikq)j^'{v + ajMkq)] [uikq)Yiv + aj5)u{kq)] 



{-g^.kl + k^k-^) {v' + a' 



(3.10) 
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Figure 5: Energy spectrum of quark-antiquark pairs emitted from the 4-body decay tr — raqq 
with an invariant mass mqq, normahsed to the 2-body decay rate, r(TR — ra). The shown quantity 
is independent of eq. The other parameters are set to nig — 1.1 my, fa — 10^^ GeV and y — I- 

where v and a are the respective vector and axial couplings and ky = kg + kg. The 
remaining integrations over cos^y and mrv, 

dr(?R raqq) ^ 1 A^/^(m|, m|, m^^^) mly - mg^- rugg 

dm„g ~ 128(27r)6m^ 2m^ 2m^y 2 

dcos6'y / dm^y 

-1 J raqq 

need to be performed numerically. Here A(a, 6, c) = + 6^ + — 2(a6 + ac + 6c), Q'y is 
the polar angle of fiy, and |7V4p is the squared matrix element with the external quark 
current replaced by expression (3.10). 

In Fig. 5 we show the obtained energy spectrum of the quark-antiquark pair stemming 
from the 4-body decay for several choices of the stau and axino masses, = l.lm^ 
and fa = 10^^ GeV. It is normalised to the 2-body decay width given by (3.1). As 
evident from the figure the effect of the Z boson resonance becomes important only for 

— rria ~ 100 GeV. 

Since pions decay fast, the major impact of hadronic energy release on primordial nu- 
cleosynthesis stems from nucleons that hadronise from the qq pair (see Sec. 4.3). Therefore 
we consider only qq pairs with nigg > 2 GeV. Moreover, cc and bb pairs only contribute 
for invariant masses above the J/^(1S) resonance with rrij/^ is — 3-0 GeV and the T(1S) 
resonance with mx,is = 9.5 GeV, respectively. The contribution from ti pairs is strongly 
suppressed due to the large threshold m^i ^ 350 GeV. In the matrix element all quarks are 
taken to be massless. This approximation is justified by the behaviour of the cross section 
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Figure 6: The discrepancy between the LL approxmiation (3.12) and full expression (3.1) for the 
2-body decay of the stau NLSP, 1 — r(TR Ta)LL/r('?k ^ ™)j for mj^^ = 1.1 rur (solid lines) and 
TOjj = 1.01 TO?; (dotted lines, red). We have set |eQ| — 1/3, y = 1 and taken the limit rn^/m- <C 1. 



for e^e — ?■ gg as a function of the centre of mass energy, which is well described by step 
function thresholds for the heavy quarks. 

3.4 Lifetime of Charged NLSP Sleptons in Axino Dark Matter Scenarios 

The results on the branching ratios of the 3- and 4-body decays obtained in Sees. 3.2 
and 3.3 show that the stau NLSP lifetime = l/L^ot t)e very well estimated by the 
2-body partial decay width alone, given by (3.1): r(rR — ra). Indeed, this decay 

width is dominated by the leading logarithmic (LL) term of the 2-body decay, 

81a^ e% m- on o f 2 

The mass of the axino enters only through a kinematic factor in this leading term, and the 
lifetime quickly becomes independent of if the latter is sufficiently smaller than m^f, 
i.e. when Aa is small enough. Moreover, it is useful to look at the LL expression (3.12) to 
see analytically the dependency of the lifetime on the parameters of the underlying model. 
It is highly sensitive to the charge of the heavy (s)quarks, oc I/cq, and less so to the 
coupling y. 

Figure 6 shows the discrepancy between the LL approximation (3.12) and the full 
expression (3.1), 1 — r(rR — t- Ta)i^i^ /V (jyi — t- to). This is maximised for small fa and m^, 
and large and but small rrij^/m^. We can also compare our LL result (3.12) to 
Eq. (2) in Ref. [37], where a cut-off method [32] had been used to calculate the 2-body 
decay. Fixing the mass scale to m = 100 GeV as suggested in [37], the discrepancy is 



-13- 




mf[TeV] 

Figure 7: The lifetime of the stau NLSP as a function of nir for m~/m~ ^ 1, jegj = 1/3, y = 1 
and fa values ranging from 10^° to 10^"' GeV. The solid, dashed and dotted lines indicate for 
nig — 1.1 m^:, rrij^ — 2 my and tojj = 1.01 m^^, respectively. 

insensitive to fa as, in the leading logarithm, it enters in the same way in both cases. 
However, due to this fixed mass scale used in [37], its discrepancy to the exact full result 
(3.1) can be much larger than that of the LL result (3.12) in our calculation. For this 
case, the error increases with increasing and m^, and decreasing /„. In general, it 
represents an overestimation of the exact result and the discrepancy reaches more than 
35% for m:^: = 1 TeV, rn§ = 2mf and fa = 10^ GeV. 

The stau lifetime (including the non-LL terms) is plotted in Fig. 7 as a function of 
its mass m:^ . It illustrates the sensitivity of the lifetime not only to nir but also to fa 
and m^. We have taken the limit of m-fm- <^ 1 where the effect of a changing 
becomes negligible. For stau masses in the range that would be accessible at the LHC, and 
the allowed fa range (1.1), the lifetime can range from about 10~^ s to 10^ s. A lifetime 
larger than 10^ s can be excluded by catalysed big bang nucleosynthesis (CBBN) as will 
be discussed in Sec. 4.3. In the analyses that follow, we use the full expression of the 
2-body decay to calculate the stau NLSP lifetime. For demonstrative purposes we often 
set I eg I = 1/3 and y = 1. We also only show our results for > 80 GeV [71], the limit 
set from searches for long-lived staus at the Large Electron Positron (LEP) Collider. 

4. Cosmological Constraints on Axino Dark Matter Scenarios 

In this section we use the results for the stau NLSP decays obtained above to systematically 
study cosmological constraints on axino dark matter scenarios. First we focus on the relic 
axino density, which is constrained by the present dark matter density. Then we will explore 
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constraints from structure formation, where upper limits on the present free-streaming 
velocities of axino dark matter emerge in order to respect the power spectrum observed 
at small scales. Lastly, we consider BBN which can considerably constrain the parameter 
space of the considered SUSY hadronic axion models. We will close this section with a 
comparative discussion of the cosmological constraints. 

For concrctcncss, we focus also in this section on the stau NLSP case. Nevertheless, 
all addressed constraints - except those associated with photodissociation, which require a 
minor modification but will turn out to be subleading anyhow - apply also directly to the 
alternative eR and /xr NLSP cases after obvious substitutions such as 171^ i^en/jlR- 

For the following considerations, we assume a standard thermal history with a post- 
inflationary reheating temperature in the range of T[ < Tr < fa, where T[ ~ is the 
freeze-out temperature of the staus, that is, the temperature at which the staus decouple 
from the thermal bath in the early universe. By considering Tr < fa, wc focus on scenarios 
in which the PQ symmetry is not restored during or after reheating. In fact, we assume 
that the PQ symmetry was broken before inflation and not restored afterwards. For large 
fa, the relic axion density is then governed by the initial misalignment angle G, of the 
axion field with respect to the CP-conserving position; see [6,7] and references therein. This 
allows us to keep the presented dark matter density constraints conservative by assuming 
fla S^dm which is possible even for fa above 10^^ GcV since Gj can be sufficiently small. 
We also assume implicitly properties of the saxion - the bosonic partner of the axino that 
appears in addition to the axion - such that its potential cosmological effects [16-19,72] are 
negligible in this study. In more general settings, the saxion can be a late decaying particle 
and as such be associated with significant entropy production [16-19] and thereby with a 
non-standard thermal history. This can affect the cosmological constraints [72] including 
those derived below, which rely on the hypothesis of a standard thermal history. 

4.1 Constraints from the Present Dark Matter Density 

In this section we discuss constraints imposed by the dark matter density [71,73] 

nl^h' = 0.105t'ofsl, (4.1) 

as obtained from the measurements of the cosmic microwave background (CMB) anisotropics 
by the Wilkinson Microwave Anisotropy Probe (WMAP) satellite, with h = 0.731qq3 de- 
noting the Hubble constant in units of lOOkms"^ Mpc~^. Note that this nominal 3a range 
is derived assuming a six-parameter flat ACDM model. If the fit is performed in the context 
of a more general model, a larger range is possible even with additional data from other 
cosmological probes [74].^ The dark matter density (4.1) constrains the relic axino density, 

n~a = f)f + n^J^ < f)d.n - na, (4.2) 

where we include contributions of thermal relic (therm) or thermally produced (TP) axinos 
and of non-thermally produced (NTP) axinos from the stau NLSP decays discussed above. 

^The updated values extracted from WMAP 5- and 7-year data [75, 76] have also been published and 
are consistent with the above range inferred from the WMAP 3-year data set [73]. 
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The freeze-out temperature of the axino can be very high due to its extremely weak 
interactions. For example, ~ 10^ GeV (10^^ GeV) has been found for /„ = 10^^ GeV 
(10^^ GeV) [22,33,77], which is thereby similar to the one calcuated for axions [6,8,78]. 

For a reheating temperature Tr > T", axinos were initially in thermal equilibrium with 
the primordial plasma, and decoupled while still relativistic, at S> m^. Accordingly, 
the resulting thermal relic axino density is given by [22,31,33,77] 




« ma/{2 keV) , (4.4) 

where g^siTf") is the number of effectively massless degrees of freedom at decoupling. As 
is typically very high, we can consider all particles in the MSSM as relativistic at this 
time, and with the addition of the axion and the axino, we can use 5*5 (Tf) « 230, which 
leads to (4.4) above. The constraint (4.2) thus translates into < 0.2 keV [22,31] even 
in the most conservative case with ^-'^^ + ^ ^^dm- As we will see in Sec. 4.2, this 
places thermal relic axinos in the 'hot' dark matter category, where even more restrictive 
constraints apply. It should however be noted that an axion cold dark matter (CDM) 
scenario, ~ ^dim in which thermal relic axinos with <^ 0.2 keV provide a tolerable 
contribution to hot dark matter is a viable option [33]; see also Sec. 4.2 below. Interestingly, 
most of the BBN constraints discussed in Sec. 4.3 will also apply to such scenarios with a 
very light axino LSP and a thermal relic stau NLSP. 

For Tr < T^, axinos were never in thermal equilibrium with the thermal bath but 
can still be produced via scattering of coloured particles in the hot MSSM plasma [31, 
33,36,77,79]. Working within the framework of SUSY QCD, hard thermal loop (HTL) 
resummation and the Braaten-Yuan prescription [80] have been used to calculate the rate 
of this thermal axino production in a gauge-invariant way and consistent to leading order 
in the strong coupling as = 5's/47r [33]. Considering Tr as the initial temperature of 
the radiation-dominated epoch, this rate leads to the following expression for the relic 
density [33,35] 

.|^.-...j™.(i^)(^)(^)^(^), ... 

where we evaluate gg (Tr) = ^/47^a^iTR) according to its 1-loop renormalisation group run- 
ning within the MSSM from as(mz) = 0.1176 at mz = 91.1876 GeV. We work with (4.5) 
since it emerges from a gauge-invariant treatment that is consistent to leading order in gg. 
Note however that this treatment relies on the weak coupling limit, gg <^ 1, and thus on 
high temperatures, T ^ 10^ GeV [33].^ 

By confronting fl~^ with the constraint (4.2), one obtains upper limits on Tr for given 
rria and /„ (or equivalently lower limits on /„ for given and Tr) [31,33-36,72,81]. 

*For alternative approaches that consider also thermal axion production at lower temperatures, see [36, 
79] where [36] also includes an additional squark-squark-gluino-ax;ino vertex and a difTerent Tr definition. 
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For (4.5) and Jlf^P ^ ^ Q^^^^ 

we have shown these hmits in Fig. 1 of Ref. [35]. For 
example, insisting on a high reheating temperature of Tr > 10^ GeV, which is desirable 
when considering the viability of standard thermal leptogenesis, and on axinos being cold 
dark matter, > 100 keV (explored more in Sec. 4.2), we are forced to consider fairly high 
fa values of /a > 3 X 10^^ GeV [35]. Note however, that in this paper we do not constrain 
ourselves to any specific baryogenesis scenario or range of fa values other than (1.1). 

Along with the thermal relic or thermally produced axinos, there are also non-thermally 
produced axinos from out-of-equilibrium NLSP decays. Due to the highly suppressed 
coupling of the axino to MSSM particles, the stau NLSP typically decays after its freeze- 
out, i.e., when T <^ T[. We thus have to consider only the yield of the stau NLSP after 
freeze-out, Yf = nf/s, where s is the total entropy density of the Universe and nf the 
number density that the stau NLSP would have today, if it had not decayed. As each stau 
decays into one axino, the relic density of the axino from this non-thermal production is 
then given by the following expression [28-32] 

n^'^^h'' = maYrs{To)hyp,, (4.6) 

where Pc/[s{TQ)h'^] = 3.6 eV as obtained from the critical density Pc/h!^ = 8.1 x 10"^''' GeV^, 
the present temperature Tq = 2.73 K = 0.235 meV and the number of effectively massless 
degrees of freedom governing the entropy density today g*s{To) = 3.91. 

The stau NLSP yield If can be calculated numerically and is model-dependent [82-87]. 
For the purposes of our work, with the focus on the tr NLSP setting, we consider the 
following three characteristic approximations: 

~ X 10-12 (^^) , = 0.7, 1.4, 2.0 (4.7) 

where k quantifies representative differences in the yield due to possible mass degeneracies 
of the stau NLSP to other sparticles and hence the extent of coannihilation. The value 
K = 0.7 corresponds to the case with = 1.1 and <C Tn^^fi, in which primordial 
stau annihilation involves only (anti-)staus in the initial state [82].^ The yield associated 
with K = 1.4 is encountered if there is either additional stau-slepton coannihilation corre- 
sponding to m^r < me,fj, < 1.1 [82] or additional stau-bino coannihilation corresponding 
to m^r < < 1.1 m^r (see contours close to the dashed line in the right panel of Fig. 3 
in [84]). For an approximate degeneracy of with both rrie^fi and m^, simultaneous 
stau-slepton-bino coannihilation can lead to an even larger yield with k = 2.0 in (4.7) (see 
Y^ contours close to the dashed line in the left panel of Fig. 3 of [84]). 

Putting the above together, we can rewrite (4.2) to obtain the dark matter constraint 
on the stau abundance prior to decay, IV < ^f^dm ^^^^ 

.»,.....o-"(^M^)(i2i^). 

Figure 8 shows the corresponding constraints on and that emerge when (4.7) is 



^The bino mass rtig = 1.1 my considered in Ref. [82] represents a typical mass splitting in regions with 
rtig > tUt encountered in scenarios such as the constrained MSSM (CMSSM). 
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Figure 8: Constraint imposed by h'^ + n^^^ < 0.126 for /„ = 10" GeV, Tr = 10^ 10^, 
10^, lO^GeV and Tr such that f2~^ is neghgible. The solid, dashed and dotted hues are obtained 
for the thermal relic stau yield (4.7) with k = 0.7, 1.4 and 2.0, respectively. For given Tr and k, 
the region to the right of the corresponding line is disfavoured by the dark matter constraint. 

confronted with (4.8) for n^^^'h^ = 0.126, < J^dm, fa = GeV and Tr = 10^ 
10^, 10^, 10^ GeV and a sufficiently small Tr value, where can be neglected. The 
solid, dashed and dotted lines are obtained respectively for k = 0.7, 1.4 and 2.0. The 
disfavoured region is the one to the right of the respective curve. Differences with respect 
to K become visible/significant only for mjj > 10 GeV, which illustrates that ^^-"""^ is 
negligible for 112^ ^ 10 GeV in the considered ni:^ range. By increasing (decreasing) fa 
by one order of magnitude, the Tr value at a given curve will increase (decrease) by two 
orders of magnitude. Note that the curves can also be read as the upper limits on Tr that 
are imposed by the dark matter constraint (4.2) for given m^, m^, and k when fta ^ ^^dm- 

In the coming sections, wc will take the most conservative point of view with + 
fla ^dm (i-G-, with Sufficiently small Tr and Qi) when showing the dark matter con- 
straint. By confronting (4.8) with 1~, as described by (4.7), the associated constraint 
will then disfavour and combinations that are located to the right/above of the 
grey bands in Figs. 9 and 13 below. Note however that the behaviour of the dark matter 
constraint for sizeable Tr leading to sizeable can easily be inferred from Fig. 8. 

One can also consider a more phenomenologically-drivcn approach where one assumes 
that all of the observed dark matter density is composed of axinos from NLSP decays, 
^dm — ^^a"""^- Irrespective of the freeze-out behaviour of the stau NLSP, this requires the 
following yield 

y, = 3.8xlO-f^^VH^y (4.9) 
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(and again Tr and Qi sufficiently small such that + Qa ^ ^^dm)- For small values of m^, 
typical values of the thermal relic stau yield, as described by (4.7), cannot reach (4.9). The 
higher yield can however be accounted for by additional non-standard sources of staus such 
as late-decaying scalar fields. In the coming sections, we explore cosmological constraints 
also for scenarios with 1% described by (4.9). 

Before proceeding, let us stress that there exists a lower limit on below which (4.9) is 
not consistent with the assumption of a standard thermal history. Indeed, that assumption 
implies that the energy density of the non-relativistic stau NLSPs prior to decay should 
not exceed the one of radiation: pf < ff prad with a fraction /f ~ 0.1 or smaller. Imposing 
this requirement at the time of the stau NLSP decay, t = t^, we find 

33.2 MeV / 50 \ /l/3\ V fa 



ff Vlog[?/VaV(2m|)] jKeq) \ 10^3 GeV 



( 



loT 



If y/2/ 100GeV \ / g.5(r^) \ / 3.36 
GeY) \ )\ 3.91 ) \g.{Tf)) ' ^^"'^^ 



as obtained for mg; <C nif when using the yield (4.9) and the LL approximation of Tf in the 
sudden decay approximation in which all tr's are considered to decay at Tf. Here g*{Tf) 
and 5*5 (t;?) are the effective number of relativistic degrees of freedom governing the energy 
density and entropy density, respectively, at the time the stau NLSP decays into the axino 
LSP. For example, if the staus decay after electrons and positrons become non-relativistic, 
9*s{tt) = 9*s{To) = 3.91 and g*{Tf) = 5*(ro) = 3.36. Note that (4.10) with ff ~ 0.1 is 
respected for the range considered in Figs. 9 and 14 below. 

4.2 Constraints from Structure Formation 

Let us now address the constraints on axino dark matter scenarios imposed by the matter 
power spectrum and related studies. Depending on details of their primordial origin, axinos 
can fall into the categories of cold, warm, and hot dark matter. As warm dark matter 
(WDM) or hot dark matter (HDM) , they can be associated with significant suppression of 
structures on scales below a potentially sizeable comoving free-streaming scale Aps • Studies 
of cosmic structures can thus provide upper limits on Aps or, equivalently, on the present 
free-streaming velocity Vpa^ of axino dark matter. We compile various limits on fpg obtained 
in WDM studies and confront those limits with the free-streaming velocities obtained for 
thermal relic axinos, thermally produced axinos, and axinos from stau NLSP decays. 

In WDM investigations of cosmic structure formation, it is typically assumed that all 
dark matter consists of one species with mass mwDM that was once in thermal equilibrium 
with the primordial plasma and freezes out while relativistic at T^^^ ^ w-WDM- Corre- 
spondingly, for a Majorana fermion of spin 1/2, $7-wdm/''^ is given by (4.3) after the obvious 
substitutions, where g*s{T^^^) is now fixed by the requirement riwDM = ^dm for a given 
"^WDM- With an initial root mean squared momentum at freeze-out of {p{T^^^)) = 
3.151 TfWDM^ ti^g present root mean squared free-streaming velocity is governed by mwDM^ 

. rms,O^WDM _ n 75 ^ ^Wi^\ ( 100eV\ , . 

^ s [0.105) UwDMy' • ^ ' 
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Probe 


™min 
'"'WDM 

[keV] 


/ -7 1 ' 1 rnax 
V'^FS 

[km/s] 


(Afs)"'^" 
[Mpc] 


Reference 


Number of MW satellites 


2.1 


0.01 


0.2 


[91] 


Luminosity function of MW satellites 


1 


0.03 


0.4 


[92] 


Phase space density in dSphs 


0.7 


0.06 


0.5 


[93] 


Phase space density in dSphs 


0.6-1.5 


0.02-0.07 


0.2-0.7 


[94] 


Lyman-a forest {z ~ 3) 


0.75 


0.05 


0.5 


[95] 


Lyman-a forest (z ~ 2.5) 


0.55 


0.08 


0.7 


[96] 


Lyman-a forest (2 < z < 4) 


2-2.4 


0.01 


0.1 


[97, 98] 


Lyman-a forest (2.0 < z < 6.4) 


4 


0.005 


0.07 


[99] 



Table 2: A selection of constraints on WDM particles from observations and numerical simulations. 

The lower limits on mwDM are taken from the corresponding references. The values for (^^'p™^' )"^^'^ 
and (Aps)™'*'^ are derived from m'^^^^ using (4.11) and (4.12), respectively, with Zcq ~ 3200 and 
teq — 1.8 X 10^^ s. The limits from [91,96-99] were quoted as 2a hmits, and a 3a difference with 
respect to observations was reported in [95], while confidence levels were not quantified in [92-94]. 

Using this expression, we obtain the upper limits on Upg listed in Table 2 from the respective 
WDM constraints provided in the literature in terms of a lower limit on mwDM- In Table 2 
we give each limit also in terms of an upper limit on the free-streaming scale 

with the cosmic scale factor a(t), the dark matter velocity v{t) = \p{t) \ / ^ \'p{t) \ ^ + ?n-\vDM 
and its present free-streaming velocity t;pg, which is assumed to be non-relativistic. The 
time and redshift at matter-radiation equality are given respectively by Zeq 3200 and 
feq — 1-8 X 10^^ s for a matter density parameter of Q,-aih? = 0.133 [10]. Note that the 
relation (4.12) between Aps and Upg is obtained in the limit that the considered dark 
matter population was produced (or decoupled from the primordial plasma) at ti much 
before radiation-matter equality, i.e., for ti <^ t^q [88-90]. 

Table 2 lists representative limits inferred from various astrophysical observations and 
studies of structure formation in WDM cosmologies.^ Included are constraints from the 
number of Milky Way (MW) satellites [91] and their luminosity function [92]. Here insights 
into the number and properties of satellite galaxies around our MW are taken into account 
which became accessible due to the results of the Sloan Digital Sky Survey (SDSS) only 
recently. Another source for the Ji^^^M constraint is the maximum observed phase space 
density in dwarf spheroidal galaxies (dSphs) [93,94]. Moreover, one of the most sensitive 
probes is the Lyman-a forest. Associated studies give some of the most constraining limits 
on the thermal WDM mass [95-99] . The data are available over a large range of redshifts z 

®For a recent review with references to further WDM studies, we refer to [100]. 
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and allow us to study the power spectrum down to very small scales. In particular, Ref. [99] 
uses High Resolution Echelle Spectrometer (HIRES) data in combination with data from 
SDSS to place a lower limit on the WDM mass of mwDM ^ 4 keV. As can be seen 
in Table 2, this 4 keV limit is about a factor of two (or more) larger than the limits coming 
from other observations. The spread in the mwDM limits reflects the fact that WDM 
constraints are subject to ongoing research and to potential systematic uncertainties. In 
the following we will put some focus on 

'TZwDM ^ 1 keV = fp^^'^ ^ 0.03 km/s (conservative), (4-13) 
rnwDM>4keV = Wp^^'° < 0.005 km/s (most restrictive). (4.14) 

These two limits will also indicate the range into which most of the limits fall that are 
given in Table 2. Moreover, future studies such as investigations of WDM constraints from 
cosmic shear [101] are expected to be sensitive to mwDM limits in that region. 

For thermal relic axinos which freeze-out with a mean squared momentum of {p{T^)) = 
3.151 T", the root mean squared value of their present free-streaming velocity is given by 

/ \ 1/'^ / 

, rms,Oxtherm ^m / 230 \ / 100 eV 



Thus, the axino mass limit derived from the dark matter density constraint in Sec. 4.1, 
rria ^ 0.2 keV [22,31], implies (^'pg^''^)|'^*^™ ^ 0.3 km/s which is an order of magnitude above 
the conservative limit (4.13), i.e., thermal relic axinos would suppress small-scale structure 
too much when providing all the dark matter. They are thus considered too 'hot' and 
are disfavoured as the dominant component of today's dark matter density. Nevertheless, 
they can provide hot dark matter in addition to neutrinos and can coexist with some other 
species providing CDM such as axions [33] . For such a setting, we find 

ma < 37 eV (4.16) 

by translating the gravitino mass limit < 16 eV obtained at the 2a level for mixed 
models with CDM and thermal relic gravitinos with 5*5 (Tj^) ~ 100 in [96] to the mixed 
case with CDM and thermal relic axinos with g^,s{Tf) — 230. In fact, a comparison of the 
Lyman-a limits listed Table 2 indicates that a more restrictive limit could emerge once 
such a mixed scenario is confronted with the data sets considered in [97-99] . 

For thermally produced axions, we estimate their mean squared momentum at the 
reheating temperature Tr, where their production is most efficient, as {p{T^)) = 3.151 Tr. 
Although never in thermal equilibrium, this choice is motivated by the fact that these axinos 
are produced in scattering of particles that are in thermal equilibrium with the primordial 
plasma. We estimate {v^^' )~^ accordingly by (4.15) after the substitution T" Tr. 
As mentioned in the previous section, we can use g^s ~ 230 for large Tr. Clearly, for 
nia ^ 100 kcV, the present free-streaming velocity of these thermally produced axinos is 
well below the constraints shown in Table 2, which places them into the category of CDM. 
For 12 keV < mj < 100 keV, thermally produced axions can be considered as WDM but 



- 21 - 



will still be compatible with the most restrictive limit in Table 2: (^'fs'^'*')s^ ^ 0.005 km/s. 
The more conservative limit (4.13) would even allow for as light as about 2 keV. The 
lower limit on can become relevant for models of inflation and baryogenesis in axino 
dark matter scenarios with 2± ridm^ the lighter allowed by structure formation the 
higher will be the Tr value allowed by the dark matter constraint (4.2) for a given /„; see 
also Fig. 8 above and Fig. 1 in Ref. [35]. 

For axinos produced non-thermally in stau NLSP decays, the present free-streaming 
velocity (^^Fs)a^^ inferred from the 2-body stau decay calculated in Sec. 3. Since 

the 2-body decay is the dominant decay channel, most non-thermally produced axinos have 
an initial momentum given by 

2 2 2 

, ^ , M m~ — m~ — mt 

1^(^01 = 9 " ■ (4-17) 

By considering the redshift of (4.17) due to the expansion of the universe, one then obtains 
their present free-streaming velocity [102], 

^NTP_ To ml -ml- ml ( g.s{To)V''' ( ^t^^ g*{rr)rl~^''^ 



which is valid in the usual case that they are produced during the radiation-dominated 
epoch, ti <C teq, and that they are non-relativistic today, |ps:(to)| = ''^a('i'Fs)s'^^- Here 
Mpi = 2.4 X 10^^ GeV is the reduced Planck mass. Note that we use the sudden decay 
approximation in which all staus are considered to decay at time U = In the following, 
we will use our full result (3.1) to calculate (fFs)5^^ ^'^^ limit m-r 0. Nevertheless, 
the main features are already described by the following estimate obtained with the LL 
approximation (3.12), 

^NTP ^ . 9 km / 50 \ / V3 \ ' ( U \ ( lGeW \ 

^ • s W[yyi/{2ml)])\eQ) VlO^^GevA rn, ) 
X 1/2 / lOOGeVN (jm_\ 



-100 GeV/ 



While a related estimate was given in [102], (4.19) illustrates explicitly the dependence on 
the logarithmic factor and on eg. 

Although axinos produced in primordial stau NLSP decays have a non-thermal dis- 
tribution, we will now confront (^'fs)^''^^ with the limits (4.13) and (4.14) inferred from 
WDM studies (cf. Table 2) assuming a thermal distribution. This approach is motivated 
by the results of [88] which show that velocity limits derived with (4.11) from mwDM limits 
can, to a reasonable extent, also be applied to the case with a monochromatic spectrum. 

Figure 9 shows contours of the present free-streaming velocity of axinos from stau 
NLSP decays (^^J^s)^"^"^. For ^f^^ = ^dm, as obtained with the stau NLSP yield (4.9), 
the regions to the left of the (purple) lines labeled with 0.03 km/s and 0.005 km/s are 



'^For treatments considering the phase space distribution beyond the sudden decay approximation, see 
Refs. [28,29,31,103-106]. 
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Figure 9: Contours of the present free-streaming velocity of axinos from stau NLSP decays, 
(^Fs)a'^^' ~ ^-ini^, fa — 10^'^GeV, leg] = 1/3 and y = 1. For the case in which fldm is 

provided by f^-"""^, the conservative hmit (0.03 km/s) and the most restrictive hmit (0.005 km/s) 
hsted in Table 2 disfavour the regions to the left of the respective (purple) lines. The stau NLSP 
yield (4.9) satisfies 17-"'"^= ftdm by construction everywhere in the shown parameter space. The 
grey band indicates the region with 17-'^^= fi^^ for the yield (4.7) with k = 0.7 or 1.4. 

disfavoured by the respective limits (4.13) and (4.14). The grey band indicates the region 
in which the axino density ^^-"""^ agrees with the dark matter density range given 
in (4.1), when using the thermal relic stau yield (4.7) with k = 0.7 or 1.4, where the latter 
accounts for potential stau-slepton coannihilation. We see that scenarios within this band 
respect (^^Fs)a"'"^ ^ ^-^^ km/s (corresponding to mwDM > 2 keV listed in Table 2) and 
that the most restrictive Lyman-a forest constraint, ^ 0.005 km/s, starts to 

disfavour some of that region only for rrir ^ 3TeV. As is evident from (4.19), the velocity 
contours are sensitive to eq and fa- An increase in the charge of the heavy (s) quarks to 
I eg I = 1 reduces the velocities by one order of magnitude while an increase in fa by an 
order of magnitude shifts the contours in the opposite way by the same amount. Thus, the 
constraints become less restrictive for higher eg and lower fa values. 

4.3 Constraints from Primordial Nucleosynthesis 

For the baryon-to-photon ratio rj (= ne/n^) inferred from studies of the CMB, standard 
BBN (SEEN) predicts the primordial abundances of deuterium, helium, and lithium in 
good overall agreement with observations [10]. When dealing with particle physics scenarios 
beyond the SM, one therefore has to check for non-standard processes that may spoil this 
success. In our case, the outcome of BBN may be affected in the three following ways: 

(1) Axinos may contribute to the radiation density at the onset of BBN prior to 
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annihilation, t 1 s, in a way analogous to the case of extra light neutrino species. 
The resulting speed-up of the Hubble expansion rate would then lead to a more 
efficient output of ^He. The inferred primordial ^He abundance thereby implies upper 
limits on the equivalent number of additional neutrinos AiVj^; cf. [107-109] for recent 
studies. (Prior to 6=*= annihilation, ANi, = N,y — 3, where the effective number of light 
neutrino species iVj^ parametrises the number of light relativistic species and their 
contribution to the total relativistic energy density, i.e., p^ad = vr^p + 7Nu/A] 
for T ~ 1 MeV.) While the resulting bounds on ANi, are still somewhat uncertain, 
we will calculate ANi, potentially contributed in the form of axinos at T ~ 1 MeV. 
Moreover, we will show how an upper limit on ANi^ translates into a constraint on 
the considered ax;ino LSP scenarios with a stau NLSP described by (4.9). 

(2) In late stau NLSP decays, energetic SM particles can be generated that could scatter 
off and/or dissociate the primordial nuclei. If injected at t > 100 s, emitted hadrons 
can reprocess the primordial light elements significantly via hadrodissociation [60,61, 
65-68,70]. At later times, t > 10^ s, also photodissociation caused by electromagnetic 
energy injection can become efficient [60,61,69,110-112]. For the considered scenarios, 
we have already published first results on the associated BBN constraints from both 
hadronic and electromagnetic energy release in a Letter [39]. This section expands 
on our previous results by exploring those constraints in more detail. 

For t < 100 s, the decay products lose their energy very efficiently and are quickly 
stopped via electromagnetic interactions. While dissociation of nuclei is not impor- 
tant at these times, additional slow- moving hadrons can take part in proton-neutron 
interconversion processes [69, 112] which could increase the abundances of ^He and 
D. Without a dedicated numerical treatment of such non-standard BBN processes 
at hand, we do not consider the associated limits in this work. As in gravitino LSP 
scenarios with a stau NLSP [113], we anyhow expect those limits to be milder than 
the ones associated with hadro/photodissociation. 

(3) If the stau NLSP is sufficiently long-lived, > 10^ s, negatively charged staus can 
form bound states with the primordial nuclei and thereby catalyse the formation of 
^Li and ^Be in the early universe [45-47]. Together with restrictive observationally 
inferred limits on the primordial abundances of ^Li and ^Be, this catalysed BBN 
(CBBN) has allowed us to constrain the parameter space of the considered models 
considerably. While CBBN constraints on axino LSP scenarios with the stau NLSP 
were first presented in our Letter [35], they will also be addressed in this section. 

Let us first turn to the BBN constraint on AN^ (= A'^^ — 3) governed by the observation- 
ally inferred primordial ^He abundance (for given ry). For example, for a ^He mass fraction 
Yp = 0.240±0.006, favoured contours (95% CL) with rj/lQ-^^ = 5.7±0.8 and AT^ = 2.4±0.8 
were obtained in Ref. [107]. However, larger Yp = 0.2565 ± 0.001 (stat) ± 0.005 (syst) [108] 
and Yp = 0.2561 it 0.011 [109] have been obtained recently, which point to larger N^, such 
as = S.GStojo (9^% CL) for rj /IQ-^^ = 6.47 and a neutron lifetime of 885 ± 0.9 s [108]. 
This exemplifies the uncertainty of the BBN bound on ANi, with, e.g., AA^^ ~ 1 being still 
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compatible with observations. With this in mind, we now show that the current restrictions 
on ANi, do not impose (significant) constraints for most of the parameter space. 

The equivalent number of additional neutrinos AN,y provided in the form of axinos 
depends on their primordial origin, the axino mass and the cosmic time at which AN,^ 
is considered. At T 2± 1 MeV, prior to the epoch of annihilation, thermal relic axinos 
have the temperature = [10.75/c/*5(rf )] ~ 0.36 MeV since they decouple when 
9*s{T^) ~ 230. A thermal axino population with ^ 0.36 MeV will then be fully 
relativistic and contribute 

T~ ( 1 7^ \ 
(AAr,)!^^-- ^ _^ ^ ( J ^ ^^-^"^ 

which applies when the respective dark matter constraint and the structure formation con- 
straint (4.16) discussed above are satisfied. Since also thermally produced axinos have 
basically a thermal spectrum [33] and originate from the high-tcmperature epoch with 
g,,Q ~ 230, the associated population can again be characterized by given above. Ac- 
cordingly, for nia <C 0.36 MeV, this will be a relativistic population that provides 

(...)F.o.o.,™,„(i^)(li!!|^)^(^). 

Indeed, this contribution cannot exceed (AA,^)-'^''™ since thermal processes cannot lead 
to an axino yield that is larger than the thermal relic yield Y^^"^"^ [31,33]. Non-thermally 

produced axinos from stau NLSP decays with r^^ < 1 s can also constitute a relativistic pop- 
ulation at T ~ 1 MeV, as discussed for the neutralino NLSP case in [31]. With the 2-body 
decay being the dominant decay channel, |^%(ii)| given by (4.17), the LL approximation of 
and the sudden decay approximation, we find the associated contribution 

. 6.3 . I0-» (i^jpj^) (f )^(l54v) 

V/VlOOGeVV \ g*{TrY'^ 22) 



100 GeV/ V ) K'^^'^^ J 9*s{TrY'^ 

This expression applies to cases in which non-thermally produced axinos become non- 
rclativistic at a temperature T^r ^ 1 MeV, where refers to temperature where the axino 
momentum satisfies Pz{Tnr) = w-s- Using again the LL and sudden decay approximations, 

,2 f2 ||'^^2\^ \ / „ \ 2 



^ 1 8 X 10-^ (^^E^lfllMA\ /egV/ lO^GeV N 
ma ' \ 30 J \l/3j \ fa J 

■ lOOGeV \ / n / 3.91 \ g.sjrrf'' ^S) 



Note that the total energy density at T ~ 1 MeV receives also contributions in the other 
cases with a non-relativistic non-thermally produced axino population, Tm- > 1 MeV, or 
a non-relativistic thermally produced one, > 0.36 MeV. However, if the dark matter 
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constraint (4.1) is respected, the resulting effect on the Hubble expansion rate and the ^He 
abundance is negligible, i.e., orders of magnitude smaller than the one of (4.20). 

Considering (AN^)^^'^™^^'^^ and (AA'^^)-"'"^ for characteristic values of the thermal relic 
stau yield (4.7), we find that these contributions cannot provide a sizeable non-standard 
ANi, ~ 1. In fact, they are clearly compatible with current BBN constraints on such 
as the ones summarised above. In turn, there are presently no associated constraints. On 
the other hand, considering (AN^)-'^^ for the stau yield (4.9), we find that an upper limit 
(AiVjy)™^^ translates into the following lower limit on the axino mass: 



rria > 23.8 keV 



30 \ /1/3\V fa 



log[yVaV(2m|)]; V eg ; VlO^GeV 

4/3 



_(AA^,^)"'=i 

/ rur \ 1/2 / lOOGcV \ ( g.s{TBBN) \ g.(r^)^/^ . . 

VlOOGeV/ V rn^ j\0.105/h'^J\ 10.75 J g^s{TrYl^ ' 

This constraint applies to cases in which the stau NLSPs decay before the temperature 
2bbn ~ 1 MeV is reached and is of similar restrictiveness as (4.10) for /f ~ 0.1. For the 
neutralino NLSP case, the constraint analogous to (4.24) can be found in Rcf. [31]. 

Here it should be mentioned that ANy is constrained not only at the time of BBN but 
also much later by observations of the CMB, galaxy clustering and the Lyman-a forest; 
see e.g. [107,114] and references therein. At those later times, ultra-light thermal relic or 
thermally produced axinos and non-thermally produced axinos with sufficiently small T-^^ 
will contribute to AA'^^. Here again the observationally inferred AA^j^ limits are such that 
these contributions can be accommodated, e.g., a 95%-credible interval of 2.2 < N^y < 5.8 
is obtained in [114]. For non-thermally produced axinos, we find compatibility with those 
limits for a thermal relic stau yield (4.7) and also for (4.9) respecting (4.10) with ff ~ 0.1. 

Let us now turn to the BBN constraints imposed by potential hadrodissociation of 
primordial nuclei via hadrons emitted in stau NLSP decays. At > 100 s, the emitted 
hadrons remain energetic and can affect the primordial abundances of D, ^He, ^He/D, ^Li 
and ^Li/'^Li [60,61,69,70,115-117]. However, in the region allowed by CBBN constraints 
for typical stau yields (4.7), < 4 x 10^ s (see later), the effect on the D abundance 
through the hadrodissociation of ^He provides the most restrictive constraint. This can be 
seen, e.g., in Figs. 38-41 of [69] and in Figs. 6-8 of [60]. Accordingly, we work with the 
TNLSP-dependent upper limits (95% CL) on 

Chad = Chad ^NLSP (4-25) 

given in Fig. 9 of Ref. [90] as obtained in Ref. [69] for observationally inferred primordial 
D abundances of 

D/H|^ean = (2.781°;^^) X 10-^ (severe) , (4.26) 

D/H|high = (3.98lg:^?) X 10-^ (conservative) . (4.27) 

Note that we arc not trying to give extra credence to the latter rather high D/H value but 
consider (4.27) simply as a limiting value that leads to robust conservative limits. 
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When using the hmits on (4.25), tnlsp = ^lsp = ^ ^-nd ehad is the average 
amount of hadronic energy release per decay. As the mesons produced from the r's emitted 
in the 2-body decay tr — > ra typicaUy decay too quickly to interact hadronically, the main 
contribution to the injection of hadronic energy comes from the 4-body decay tr — > raqq 
(considered in Sec. 3.3): 

1 rm,-m-^-mr dr(?R ^ raqq) 
ehad = -^/ drngguigg — ; (4.28) 

see also Fig. 5. Similarly, the mesons from this process also decay before they can interact 
with the background nuclei. As already mentioned in Sec. 3.3, we therefore place a cut on 
the invariant mass of the produced qq pairs, nigg > m^^* = 2 GeV, so as to consider only 
those that hadronise into nucleons. While (4.28) is used with this cut value to determine 
the hadronic BBN constraints in this work, we note that a more precise calculation of 
the constraints would require, in addition to dr(rR — t- raqq) / dnigq, a treatment of the 
fragmentation of the quarks into hadrons and of the propagation of the resulting hadron 
spectra when computing the abundances of primordial light elements; cf. [59,61,69,118]. 

Figure 10 shows the hadronic energy release ehad as calculated with (4.28) (solid lines) 
along with the lifetime contours of the stau NLSP (dotted lines) for |eQ| = 1/3, y = 1 
and four different combinations of fa and mj^/m^ (as indicated). The dependency of 
Chad on "la IS as Weak as that of r^, until we go to large enough where phase-space 
suppression sets in. As fa increases to larger values, the average hadronic energy release 
Chad decreases, although not as drastically as the increase in the lifetime. This is due to 
the fact that, after normalisation of Chad to F^^, the leading dependence of Chad on fa is 
oc l/log2 [y2/a/(2?ri|)] , rather than oc 1/ fl 

Using ehad as calculated with (4.28) and the upper limits on hadronic energy release 
^haT mentioned above, we obtain upper limits on the yield of the stau NLSP prior to decay 

= &/ehad . (4.29) 

Figures 11 and 12 show those limits as imposed by the severe limit associated with (4.26) 
and by the conservative limit associated with (4.27), respectively. The contours correspond 
to the indicated ^-^ad values which range between 2 x lO"-*^^ to 10~^ with darker shading 
representing more restrictive constraints. The considered ranges of and m^, the panels 
and the parameters are set as in Fig. 10. As seen from the first three panels of each figure, 
we find that increasing values of /„ are associated with increasingly tighter constraints 
on Y^. Moreover, a comparison of the lower panels in each figure shows that smaller 
values of m^/m^ imply tighter limits Y~^^^. This results from the enhancement of ehad 
for m^/m^ 1 described above. In each panel, one also finds common features of the 
constraints. The constraints become weaker for decreasing 80 GeV, which 

results from the decrease in ehad documented in Fig. 10, and one for 100 s, which 

is due to a relaxation of the limits [69, 90] . Furthermore, since ehad and T^f: become 

insensitive to below 10 GeV, the same holds also for Y~^^. Accordingly, the Y^ limits 
in Figs. 11 and 12 shown for mg < 10 GeV apply also to cases with a much lighter axino 
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Figure 10: Average amount of hadronic energy release per stau NLSP decay Chad (solid lines), 
as governed by the 4-body dec ay tr — Taqq, and stau NLSP lifetime contours (dashed lines) 
for I eg I = 1/3 and y = 1. Starting from the top left panel and going row by row, we consider 
(/a, m^/m?)=(10i2GeV, 1.1), (3 x lO^^ GeV, 1.1), (10" GeV, 1.1) and (lO^^GeV, 1.01). 



LSP, including even those that respect (4.16). Going into the other direction, — ?• m;jr, 
Chad becomes phase-space suppressed (cf. Fig. 10). This leads to the relaxation of the 1~ 
constraints for — )• m^f that can be seen in Figs. 11 and 12. 

Now by confronting the stau yield with the maximum yield Y~^^^ shown, e.g., in 
Figs. 11 and 12, we obtain the regions in the parameter space that are disfavoured by the 
hadronic BBN constraints associated with the D abundances (4.26) and (4.27) [39]. For 
characteristic values of the thermal relic yield (4.7) and for (4.9) satisfying f^-"'"^ = J^dm by 
construction, those regions are the ones enclosed by the short-dash-dotted (blue) lines in 
Figs. 13 and 14, respectively. As a guidance, dotted lines are shown to indicate r^^ = 10^, 
5 X 10^ and 10^ s in Fig. 13 and = 10^, 10^ and 10^ s in Fig. 14. In fact, we limit our 
consideration of the hadronic BBN constraints to parameter regions with > 10^ s. 

In Fig. 13 we have chosen the parameters such that the constraints shown in the panels 
become progressively more limiting. For the generic case of = 1.1 and a thermal relic 
yield (4.7) characterised hy k = 0.7, both the severe and conservative hadronic BBN limits 
are still absent for = 10^^ GeV (top left panel). For a larger thermal relic yield given 
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Figure 11: Upper limits on the yield of the stau NLSP prior to decay, 5^ad' imposed by the severe 
limit on late hadronic energy injection associated with (4.26). The contour shading darkens with 
stronger limits. The choice of the panels and associated parameters matches the one in Fig. 10. 

hy K = 1.4, accounting for a possible stau-slepton coannihilation, more severe constraints 
are obtained. In fact, the severe one associated with (4.26) disfavours a stau mass between 
300 and 900 GeV for < 100 GeV (top right panel). Decreasing disfavours also 
larger regions of parameter space due to the increase in ehad illustrated in the lower panels 
of Fig. 10. Accordingly, for the same yield characterised by k = 1.4 but now accounting 
for a possible stau-bino coannihilation with nig = 1.01 m:^, the conservative constraints 
associated with (4.27) appear for the first time (bottom left panel). In the extreme case 
of having simultaneous bino-stau-slcpton coannihilation with = 1.01 m^:^ and k = 2.0, 
one encounters the most restrictive case shown in that figure. There the severe D limit 
disfavours almost the entire parameter space with 10^ s < < 10^ s (bottom right panel). 

In Fig. 14 both the severe and the conservative hadronic BBN constraints already 
show up at fa = 10^^ GeV and disfavour increasingly larger areas of parameter space in 
the panels for /a = 3 x IQ-^^ and IQ-^^ GeV. This is due to the rather large stau yield (4.9) 
that must be assumed at small mjj to fulfil ^f^^^ = The 771^ dependence of the 

yield (4.9) explains also the significant variation of the hadronic BBN constraints with m^. 
In contrast, those constraints do not change for < 10 GeV when varying mj in case of 
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Figure 12: The same as Fig. 11 but for tlie conservative limit on late hadronic energy injection 
associated with (4.27). 

a thermal relic yield (4.7) as can be seen in Fig. 13. 

Electromagnetic BBN constraints have been explored in our study in close analogy to 
the hadronic ones. However, we condense the discussion since those constraints become 
relevant only for t:^ > 10^ s, i.e., in a region already disfavoured by the CBBN constraints 
discussed below. To determine the electromagnetic BBN constraints, we work with the Tr- 
dependent upper limits on ^cm = Ccm Yr imposed by the primordial abundances of D and 
by ^He/D. In particular, we consider severe D|5^ and ^He/D constraints inferred from the 
respective limits given in Fig. 42 of Ref. [69] and a conservative D^^® constraint inferred 
from the respective limit given in Fig. 6 of Ref. [111]. The electromagnetic energy release 
ecm is governed by the tau emitted in the 2-body decay with an energy Ej- = (m~ — ^7^| + 
m'^)/{2m^) in the rest frame of the tr. As each r decays into at least one which does not 
interact electromagnetically, only a fraction of Er contributes to the total electromagnetic 
energy injection [90,119]. We use the conservative estimate 

m% — m~ + mi 
= 0.3 Er = 0.3 ^ ^ ° ^ , 4.30 

to avoid that the electromagnetic BBN constraints presented are overly restrictive. Pro- 
ceeding as outlined for the hadronic BBN constraints above, we then obtain upper limits 
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Figure 13: Cosmological constraints on tlie masses of tlie axino LSP and stau NLSP for fa = 
10^^ GeV, legl = 1/3, y = 1 and given by representative thermal relic abundances (4.7). Starting 
from the top left panel and going row by row, we have (m^/my, k) = (1.1, 0.7), (1.1, 1.4), (1.01, 1.4) 
and (1.01, 2.0). Contours of r^r = 10^, 5 x 10'^ and 10^ s are shown by the dotted lines. On (above) 
the grey band, D,^^^ G n^^^ i^^^^h"^ > 0.126). Hadronic BBN constraints imposed by the D 
abundances (4.26) and (4.27) disfavour the regions enclosed by the respective short-dash-dotted 
(blue) lines. Electromagnetic BBN constraints do not appear in the shown parameter space. The 
CBBN constraints imposed by the ^Li range (4.34) and the ®Be abundance (4.35) disfavour the 
regions below the long-dash-dotted (red) and the solid lines, respectively. Not shown are the 
constraints from structure formation discussed in Sec. 4.2 which apply only within the grey bands. 
For the upper (lower) panels, it can be seen in Fig. 9 (Fig. 14, bottom right panel) that only the most 
restrictive velocity limit (4.14) disfavours parts of the grey bands for m^: > 3 TeV (m^r > 1 TeV). 



^?em ^^^^ yield 1~ prior to decay. For the settings considered in Fig. 13, we find 

that the thermal relic yields (4.7) respect the obtained limits in the shown parameter 
space such that no electromagnetic BBN constraints appear in that figure; see also Fig. 5 in 
Ref. [39]. This is different in Fig. 14 where we encounter electromagnetic BBN constraints 
in the panels with /a > 3 x 10^^ GeV. There the D-imposed constraints disfavour the 
regions enclosed by the short-dashed (blue) lines and the ^He/D-imposed ones the regions 
enclosed by the double-dash-dotted (green) lines. However, we see in Fig. 14 that these 
constraints only show up in the regions well within those already excluded by the CBBN 
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Figure 14: Cosmological constraints on tlie masses of tlie axino LSP and stau NLSP for tlie 
stau yield given by (4.9) that satisfies f^-"""^ = ^dm by construction everywhere in the 
shown parameter space. Starting from the top left panel and going row by row, we have 
{fa, mg/mr) = (1012 GeV, 1.1), (3 x lO^^ GeV, 1.1), (10" GeV, 1.1) and (10" GeV, 1.01), with 
the other parameters set to jegl = 1/3 and y = 1, as in Figs. 10-13 above. Contours of = 10^, 
10'^ and 10^ s are shown by the dotted lines. The constraints from structure formation associated 
with (4.13) and (4.14) disfavour the regions to the left of the respective solid (purple) diagonal 
lines. Here the grey band is only shown to indicate the region with """^ ~ f^dm for the ther- 
mal relic yield (4.7) with k = 0.7 or 1.4 in the first three panels and with k = 1.4 or 2.0 in the 
bottom right panel. Hadronic BBN constraints imposed by the D abundances (4.26) and (4.27) 
disfavour the regions enclosed by the respective short-dash-dotted (blue) lines. The corresponding 
D-imposed electromagnetic BBN constraints disfavour the regions enclosed by the respective short- 
dashed (blue) lines and the ^He/D-imposed ones the regions enclosed by the double-dash-dotted 
(green) lines. The CBBN constraints imposed by the ^Li range (4.34) and the ^Be abundance (4.35) 
disfavour the regions below the long-dash-dotted (red) and the solid lines, respectively. 

constraints, which will be discussed next, and are absent altogether for /a < 3 x 10^^ GeV. 

For the alternative cr and /Ir NLSP cases, the electromagnetic BBN constraints will 
look somewhat different. This is expected from the different properties of the respective 
emitted lepton; cf. Sec. 2.2 of Ref. [90]. For example, in the cr NLSP case, the electro- 
magnetic energy release is eem = -E-e = — "i| -|- mg)/(2mej^) instead of (4.30) since 
the electron does not decay. However, we refrain from a further detailed assessment since 
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the resulting electromagetic BBN constraints will still be less restrictive than the CBBN 
constraints for scenarios with the e^ NLSP or the /Ir NLSP. 

Let us now address CBBN constraints on axino LSP scenarios with a long-lived stau 
NLSP [35]. At cosmic times t > 10^ s, the presence of negatively-charged r^'s can lead 
to the formation of (^Her^) and (^Ber£) bound states once the rate of photodissocia- 
tion of those bound states drops below the Hubble expansion rate. Then, depending on 
the (^HcTj^) and (^Ber£) abundances at the relevant times, which are governed by 
and Yf, primordial ^Li and ^Be production can become efficient via the following CBBN 
reactions [45-47]^ 

(^Hctr) + D ^ ^Li + f- , (4.31) 
^He + (^He?£) ^ (®Ber£) +7 , (4.32) 
(^Be T^) + n ^- ^Be + T^; . (4.33) 

Here we work with the observationally-inferred limits on the respective primordial fractions 
of ^Li [111,121,122] and ^Be [47], 

^Li/Hlobs < 10-"-10-l^ (4.34) 
^Be/H|obs < 2.1 X 10~^^ (4.35) 

and obtain the CBBN constraints by confronting 1~ with the limits in Fig. 5 of Ref. [47]. 

In Figs. 13 and 14, the CBBN constraints associated with the range (4.34) are indicated 
by pairs of long-dash-dotted (^Li, red) lines and the ones associated with (4.35) by solid 
(^Be) lines. The regions below those lines arc disfavoured by an excess of ^Li and ^Be 
above the respective limits. As seen in Fig. 13, the CBBN constraints on scenarios with 
a thermal relic yield (4.7) follow contours of constant and become independent of 
for m~/m~ <^ 1. In fact, we find that the CBBN constraints disfavour t:^ > A x 10^ s and 
that they show only a mild dependence on 1~ in contrast to the hadronic BBN constraints. 
For the yield (4.9), the CBBN constraints can disfavour even smaller of about 10^ s 
(cf. Fig. 14) as they become more restrictive towards smaller due to the significant 
increase of Y^. 

Here one should also comment on the potential interplay between late energy injection 
and CBBN. We use the CBBN limits from Ref. [47] which have been derived for a D abun- 
dance obtained with SBBN. For an increased D abundance resulting from hadrodissociation 
of ^He, CBBN of ^Li and ^Be becomes even more efficient. For ^Li, this is evident since its 
catalysis proceeds via (4.31) and since the primordial abundance of D stays significantly 
below the one of ^He at the relevant times (even for a maximum of observationally tolerable 
hadrodissociation of ^He); cf. Fig. 2.4 in Ref. [123]. For ^Be, an increased output results 
from the final step of its catalysis (4.33), which becomes more efficient since an enhanced 
abundance of D increases the number of neutrons n at the relevant times [124]; cf. Fig. 4 in 
Ref. [47]. Moreover, the debris of hadrodissociated "^He can hit ambient ^He and thereby 

®The large ^Bc-production cross section reported and used in [46,47] has been questioned by [120], in 
which a 4-body-model study is announced as work in progress to clarify the efRciency of ®Be production. 
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fuse additional ^Li [60,69,70]. The interplay of late energy injection and CBBN will thus 
lead to constraints that can only be stronger than the ones presented in this work. Aiming 
at conservative limits, this allows us to neglect those intricacies which will have to be faced 
in future refinements of the presented constraints. 

4.4 Comparison and Discussion of the Cosmological Constraints 

Let us close this section with a brief comparative discussion of the cosmological constraints 
addressed in this work. We consider again cosmological settings with a thermal relic stau 
yield (4.7) and such with (4.9) satisfying = i^dm by construction. 

When presenting Fig. 13, negligible + Cla was assumed and thereby Tr <^ 10^ GeV 
and @i <^ 1. While the shown BBN constraints are not affected by this assumption, sig- 
nificant additional parts of the parameter space can be excluded by the dark matter con- 
straint (4.2) for Tr > 10^ GeV, as explicitly illustrated in Fig. 8. With thermally produced 
axinos dominating ridm, which are compatible with cosmic structure formation down to 
as low as 12 keV or even smaller (cf. Sec. 4.2), the structure formation constraints shown in 
Fig. 9 do not apply. In other cases with Tr <C 10^ GeV and a thermal relic stau yield (4.7), 
those constraints apply within the grey bands shown in Fig. 13. However, they arc still 
relatively mild. For the upper (lower) panels, it can be seen in Fig. 9 (Fig. 14, bottom 
right panel) that only the most restrictive velocity limit (4.14) disfavours parts of the grey 
bands for > 3TeV (m^ > 1 TeV). Moreover, the contributions to AN^ at T ~ 1 MeV 
are typically very small and clearly compatible with current limits for thermally produced 
axinos and for those from decays of thermal relic staus with 1~ described by (4.7). 

In Fig. 14, in which 1^ given by (4.9) is used, J^-"'"^ = ^dm is considered everywhere 
in the shown parameter space. This also implies negligible $7-^ + $7^ and thereby Tr <^ 
10^ GeV and ©j <C 1. Moreover, here the constraints from structure formation associated 
with (4.13) and (4.14) apply not only in the grey bands but disfavour the regions to the left 
of the respective solid (purple) diagonal lines. Accordingly, even the conservative structure 
formation constraint associated with (4.13) is significantly more restrictive than both the 
constraint (4.10) with /f ~ 0.1 imposed by our working hypothesis of a standard thermal 
history and the BBN constraint (4.24) with (AiV^)"^^^ = 0(1) at T ~ 1 MeV. In the 
regions with larger that are allowed by the respective structure formation constraint, 
the hadronic BBN constraints are often the most restrictive ones. While those become 
weaker for — > m:p, the CBBN constraints become increasingly more restrictive in that 
region which is associated with larger r:^ values due to phase-space suppression of the stau 
decays. As already mentioned, the electromagnetic BBN constraints show only up for 
> 10"^ s and only in parts of the parameter space already disfavoured by the CBBN 
constraints. We also show grey bands in Fig. 14 but merely as a guide to show where the 
stau yield (4.9) coincides with that of the stau thermal relic yield (4.7) with k = 0.7 or 1.4 
(first three panels) and with k = 1.4 or 2.0 (bottom right panel). Those are the regions 
in which one might consider (4.9) to take on its most natural values, i.e., values motivated 
by thermal-freeze-out considerations. 

Looking once more at Fig. 13, one finds that the hadronic BBN constraints are much 
more sensitive to the stau yield and the m^/m^ ratio than the CBBN constraints. While 
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the CBBN constraint disfavours basically < 500 GcV in all four panels, the hadronic 
BBN constraint may not show up at all or may disfavour even significantly larger val- 
ues of 1 TeV or above, depending on 1~ and the m^/m^ ratio. Indeed, this sensitivity 
motivated us to show those four combinations of k and m^/m^ accounting systematically 
for the possibilities of no coannihilation, stau-slepton coannihilation, stau-bino coannihila- 
tion, and simultaneous stau-bino-slepton coannihilation. While fa is fixed to 10^^ GeV in 
Fig. 13, an illustration of the fa dependence of the CBBN constraints can be found in Fig. 4 
of Ref. [35] and in Fig. 5 of Ref. [39] which also includes hadronic and electromagnetic BBN 
constraints. For [egj = 1/3 and y = I, those figures demonstrate that BBN constraints 
start to become an issue for > 10^^ GeV while settings with smaller /„ are free from 
those constraints when considering thermal relic stau yields and the limit m:^ > 80 GeV 
from searches at LEP [10]. With ehad/cm being independent of eg and oc l/cg, even 
settings with fa up to 10^^ GeV will be free from BBN constraints if jegl = 1 and y = 1. 

Finally, there is the possibility that very light thermal relic axinos may exist as HDM 
while CDM is provided by axions [33] . Here the obtained limit < 37 eV from structure 
formation constraints is significantly more restrictive than the previously known limit < 
200 cV imposed by the dark matter constraint (4.2). Interestingly, also for very light axinos 
the hadronic, electromagnetic and CBBN constraints apply since they are sensitive to the 
yield of the stau NLSP, its SM decay products and its lifetime, which is independent of 
in that regime. As mentioned, for given by (4.7), the location of those constraints can 
then be seen directly in Fig. 13 as they agree with the ones shown for < 10 GeV. 

5. Collider Phenomenology of Axino Dark Matter Scenarios 

Cosmological considerations can significantly constrain the parameter space of SUSY axion 
models. Nevertheless, in order to investigate and to determine the properties of axinos, 
one must look towards laboratory experiments. Due to the axino's highly suppressed 
interactions, its production cross sections at various colliders are very small. Instead, one 
would look for the NLSPs, the staus, either from direct pair production or from cascade 
decays of heavier sparticles. In fact, a first hint towards axino dark matter would come 
from the appearance of quasi-stable staus at the Tevatron, the LHC or future colliders.^ 

The most straightforward property of the stau to determine would be its mass m^, 
which can be measured using time-of-flight (TOF) data from the muon chambers. At 
the LHC, it is estimated that this can be done with < 1% accuracy [131, 132]. In order 
to measure other properties, such as the stau lifetime and the axino mass mg, it is 
important to be able to study the stau decays. The approach that one should take to 
analyse these decays depends on the decay length cr^r and on the size of the detector L. 
Here it is natural to distinguish the following three cases, which we now relate to the 
considered axino dark matter scenarios: 

• If the lifetime is very short such that ct^ <^ 0(cm) then the staus would effectively 
decay immediately after production. This decay length corresponds to ^ 10~^° s. 

^Collider studies of long-lived charged massive particles (CHAMPs), a category to which the staus 
belong, can be found for example in Refs. [125-130]. 
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For the SUSY axion models considered in this paper, a stau with such a short lifetime 
can only be achieved with a very low fa and a very heavy stau, as well as [egl larger 
than 1/3. For example, with the lowest fa value respecting (1.1), /„ ~ 6 x 10^ GeV, 
one would need > 2 TeV, along with m^/m^ = 2, rria/m^f <C 1, ?/ = 1 and 
\eQ\ = 1, to obtain r^^ < 2 x 10~^ s. This very high stau mass is out of the reach of 
current and near future colliders. If one considers a stau mass that is more likely to 
be probed at colliders, = 300 GeV while keeping /„ ~ 6 x 10® GeV, along with 
the selection of m^/m^ = 1.1, nia/m^ <C 1, y = 1 and IcqI = 1/3, one obtains a 
lifetime of t:^ 10""^ s and cr^ ^ 30 km. 

• If the lifetime is such that 0(cm) < ct^ < L then the staus will decay with a 
displaced vertex within the beam pipe or the detectors. With a typical detector size, 
L < 0.1 km, it is evident from the arguments given above that this case still requires 
a very low fa and a very heavy stau to be relevant for the considered axino LSP 
scenarios. 

• If the decay length is cr^ > L, there would be only few in-fiight decays within the 
detector and most of the staus would escape the detector. This is the case on which 
we focus in this section. It applies to axino LSP scenarios respecting (1.1) and with 
mf such that the stau NLSP can be probed at current or near future colliders. 

There are different ways to access stau decays in the considered third case. For r^j: ~ 
10~(^~5) s (still corresponding to decay lengths well over L), it has been argued that a 
substantial number of in-flight decays can still be observed within the detectors [133,134]. 
Also for larger r:^, decays coTild be accessible since the produced staus will be slowed down 
(due to ionisation energy loss) so that some of them are expected to be stopped within the 
main detectors [43,135,136] or in an additional dedicated stopping-detector [40, 41, 44]. 
By recording the stopping point/time with those devices, it seems possible to study the 
subsequent decays of the trapped staus and to determine their lifetime. For very long 
lifetimes, there was also a suggestion to place a water tank around the detector to trap the 
staus, such that the water can be moved to a different location at a later date for further 
analyses [42]. Even the surrounding rock has been suggested to be used for the study of 
stau decays [137]. 

In the coming parts of this section we first explore the prospects of stopping staus 
and studying their 2-body decays with the main detectors at the LHC and the ILC, i.e., 
without additional stopping-detectors. Four benchmark scenarios are proposed and studied 
for this purpose. Thereafter, we turn to the 3-body stau decays into the LSP, a tau and a 
photon, whose study might require additional stopping-detectors. Based on those decays, 
we address the possibility to distinguish between the axino LSP and the gravitino LSP at 
colliders. 

^'^Note that this method is not limited to staus with long Ufetimes. It can also be used to study staus 
with very short lifetimes but highly boosted momenta. 
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Scenario 


mo [GeV] 


mi/2 [GcV] 


tan /? 


[GeV] 


[GeV] 


m-ep [GeV] 


K 


A 


170 


560 


30 


210.3 


232.4 


270.1 


0.73 


B 


50 


550 


8 


207.6 


226.8 


212.4 


1.34 


C 


138 


435 


25 


174.3 


177.7 


215.6 


1.42 


D 


135 


700 


6 


290.2 


292.6 


292.8 


1.94 



Table 3: Benchmark points in the CMSSM that emulate some of the scenarios explored in Sec. 4. 
For all four point s, Aq = and sgn(/i) = 1. The weak scale masses of the lighter stau, the bino-like 
neutralino and the lighter sclcctron and the stau yield Yf given by the respective k with (4.7) are 
calculated with the tools SuSpect [138] and micrOMEGAs 2.4 [139,140] for nit = 172.5 GeV. 

5.1 Probing 2-Body Decays of Trapped Staus at the LHC and the ILC 

To illustrate the potential of stopping staus and analyzing their 2-body decay into the 
axino LSP at the LHC and the ILC, we consider the four SUSY benchmark points specified 
in Table 3. To simplify matters and minimise the number of parameters that have to be 
specified, we work in the framework of CMSSM. Here the gaugino masses, the scalar masses 
and the trilinear scalar interactions are assumed to take on the respective universal values 
rnx/2! "f^o and Aq at the scale of grand unification Mqx't — 2 x 10^^ GcV. Specifying also 
the mixing angle tan /? in the Higgs sector and the sign of the higgsino mass parameter 
fi, one can calculate associated mass spectra at the weak scale and stau yields prior to 
decay. For the considered benchmark points, we obtain those quantities by employing the 
spectrum calculator SuSpect [138] and MicrOMEGAs 2.4 [139,140] for mt = 172.5 GeV. 
Those tools are also used to verify that the accelerator constraints are respected for our 
chosen points. 

We have chosen the benchmark points such that they fall into the classes for which 
we studied cosmological constraints in Sec. 4. This can be seen in Table 3, where k values 
specify the associated stau yield given by (4.7) with the respective values. Scenario A 
shows no coannihilation, B stau— slepton coannihilation, C stau-bino coannihilation and D 
simultaneous stau-slepton-bino coannihilation; cf. text below (4.7). These four scenarios 
come with a bino-like lightest neutralino and with lighter staus that have < 10% left-handed 
component. This left-admixture can be reduced further by fine-tuning the parameters or 
considering other SUSY scenarios that offer greater flexibility in choosing parameters. In 
fact, we use our result for the stau lifetime derived for right-handed staus in Sec. 3 also for 
the four benchmark points even if they show a non-negligible left-right mixing. 

The SUSY hadronic axion model is not part of the CMSSM and the associated pa- 
rameters have to be specified in addition to those in Table 3. For concreteness, we choose 
mz = 10 GeV, fa = 10^^ GeV, y = 1 and [cq] = 1/3. While this scenario respects all 
constraints from structure formation and BBN, the dark matter constraint limits the max- 
imum reheating temperature to a relatively low value of about lO"' GeV; cf. Fig. 8 and 
associated explanations given in Sec. 4.1. For smaller and/or larger fa, this Tr limit 
would be more relaxed. Note also that a change in will not noticeably affect the stau 
lifetime and the associated error derived below, as long as rn~/m~ <C 1. However, the same 
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Figure 15: Expected |pf j/m^r distributions of lighter staus at the LHC with y/s — 14 TeV and an 
integrated luminosity of £ = 10 fb^^ for the CMSSM points A-D specified in Table 3. The events 
have been generated with Pythia 6.4 [141] and the cuts 1-4 outlined in [135] have been applied. 



is not true for a change in /„, on which the stau lifetime strongly depends; see (3.1). 

Let us now turn to the LHC. As mentioned above, we concentrate on the fraction of 
staus that are produced with small enough values of = |p~|/m^; such that they lose 
enough kinetic energy through ionisation to be stopped within the main detectors. For 
this setting, the authors of [135] suggest the possibilities of having a trigger for a beam 
dump should a stau be detected to have stopped, such that its subsequent decay can be 
observed, or of studying the decays of the staus during beam downtime. 

Our LHC study follows closely the approach of Ref. [135]. Using Pythia 6.4 [141] 
and the cuts 1-4 outlined in [135], we have generated proton-proton collision events with 
a centre of mass energy ^/s = 14 TeV and an integrated luminosity £ = 10 fb^^ for 
each CMSSM point given in Table 3. The resulting ]p~]/m^^ spectra are shown in Fig. 15. 
Following the strategy outlined in [135], we consider only the staus produced with (3^ < 0.45 
as those that get trapped in the calorimeters of the ATLAS detector. For the considered 
benchmark scenarios, we list the corresponding numbers of staus expected to be stopped 



^^Another recent LHC proposal [136] suggested to measure the lifetime of staus trapped in the ATLAS 
detector during scheduled shutdown periods. It is thus more suitable for studying staus with a very long 
lifetime of several days or more. However, in the scenarios that we are considering, a lifetime longer than 
about three hours is disfavoured by the CBBN constraints discussed in Sec. 4.3. 
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Scenario 


"^UR [GeV] 


m-g [GeV] 


O-SUSY [pb] 


^stopped 

T 


Tr[s] 


A 


1127.2 


1272.8 


0.66 


119 


5.4 ±0.8 


B 


1097.5 


1250.3 


0.76 


182 


5.8 ±0.6 


C 


896.4 


1008.3 


2.60 


522 


10.8 ±0.7 


D 


1372.1 


1563.6 


0.185 


58 


2.5 ±0.5 



Table 4: Expected number ]Sf^^°PP^'^ of lighter staus to be stopped in the calorimeters of the 
ATLAS detector at the LHC with ^/s = 14 TeV and an integrated luminosity of £ = 10 fb~ 
for the CMSSM points A-D specified in Table 3, and the resulting expected statistical 68% CL 
accuracy of the determination of for = 10 GeV, fa = 10^^ GeV, y ~ I and jegl = 1/3. For 
each CMSSM point, we also list the weak scale values of rriu^, rrig and ctsusy, as obtained with the 
spectrum calculator SuSpect [138] for mt = 172.5 GeV and with Pythia 6.4 [141]. 

in the ATLAS calorimeters, JV-*"^^^*^, together with the right-handed up squark mass m^^, 
the gluino mass rrig and the total SUSY production cross section (Tsusy in Table 4. Clearly 
jystoppcd Yarger for larger ctsusYj^^ and thus governed by rriu^ and rrig at the LHC. 

The lifetime of the trapped staus is determined by measuring the time between when 
the stau is stopped and when its subsequent decay is detected, where we refer to [135] 
for details on the different techniques one would use for different r^f. By extrapolating the 
results of [135] to our parameter points, we obtain the expected statistical 68% CL precision 
of the Tr determination listed in Table 4. Clearly the statistical accuracy depends on the 
number of events that can be observed, and at the LHC, this is very much dictated by the 
mass spectrum of the SUSY scenario that is realised. With the conditions assumed for the 
generation of the above events, light squarks and gluinos with masses rriq^g < 1 TeV would 
be preferred for a statistical 68% CL accuracy of the r^r determination of < 10%. 

A similar analysis can be done for a future linear collider. In Rcf. [43] it has been 
proposed to use the Large Detector Concept (LDC) [142] at the ILC to study the decays of 
the trapped staus. The staus produced in the collisions can become trapped either in the 
hadronic calorimeter or the yoke of the LDC. We follow very closely the strategy of [43] to 
determine the statistical accuracy with which the properties of the stau and the axino can 
be determined. 

Again using Pythia 6.4 [141], we now generate electron-positron collision events with 
a centre of mass energy of ^/s = 600 GeV and an integrated luminosity of £ = 250 fb^^ for 
each CMSSM point given in Table 3 . The obtained |p~|/m^ distributions for the produced 
lighter staus are shown in Fig. 16. In addition to the total distribution (solid line, red), 
we illustrate the contributions of the following processes to the production of lighter staus: 
e+e" ^ eRe^ (plus signs, green), /xr/I^ (crosses, blue), tit* (stars, magenta), XiXi (boxes, 
cyan), eRe£ and e^^CL (circles with crosses, yellow), Tir| and r*r2 (open circles, black) and 
X1X2 (filled circles, orange). This breakdown of the contributing processes allows us to see 
effects associated with small mass differences between the lighter stau r]; and the sleptons 

^•^Table 4 lists leading order results for (Jsusy- The values increase by factors of 1.2-1.4 when including 
the next-to-leading-order SUSY-QCD corrections so that j\[^^°pp'^'^ may be above the given estimates. 
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Figure 16: Expected Ipkl/mr distributions of lighter staus (solid line, red) in electron-positron 
collisions at the ILC with y/s = 600 GeV and an integrated luminosity of £ = 250 fb~^ for the 
CMSSM points A-D specified in Table 3. The events have been generated with Pythia 6.4 [141]. 
There are contributions from e+e" — >• e^e"^ (plus signs, green), JirJi^ (crosses, blue), rirj* (stars, 
magenta), XiXi (boxes, cyan), eRe£ and eJj^eL (circles with crosses, yellow), titJ and t*T2 (open 
circles, black) and X1X2 (filled circles, orange). 



' /^R and/or the lightest neutralino Xi^ Pronounced peaks just below the ones from 
direct r}*"* production come from cascade decays of e[^^ and Jl^^ at point B (top right 
panel), of x? ^it point C (bottom left panel) and of e[^\ Jl^^ and Xi point D (bottom 
right panel). In addition, one finds also for points A and B that a significant fraction 
of the staus is produced in decays of the lightest neutralino Xi- This indicates that a 
lightest neutralino with a mass close to the one of the lighter stau would be beneficial for 
phenomenological studies. 

Following Ref. [43], we have calculated the numbers of staus expected to be stopped 
in the hadronic calorimeter, N-^'^^, and in the yoke, N~°^'^ , of the LDC. For the CMSSM 
points A-D, the results are tabulated together with the total SUSY production cross section 
CSUSY in Table 5. The respective fSj ranges for which the staus are expected to get stopped 
in those detector parts depend on and are typically 0.3-0.4 (heal) and 0.4-0.5 (yoke) 
for our benchmark points; see also Table 2 in Ref. [43]. 

Note that iV-^*^' and N-"^'^ depend not only on ususY but are also very sensitive to 
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Scenario 


0"SUSY [fb] 


T 


^yokc 

T 


nir [GcV] 


Tr [s] 


ms [GeV] 


A 


84 


508 


1362 


210.3 ±0.2 


5.4 ±0.2 


10 ±64 


B 


190 


681 


1311 


207.6 ±0.2 


5.8 ±0.2 


10 ±58 


C 


310 


1569 


2933 


174.3 ±0.1 


10.8 ±0.3 


10 ±35 


D 


3.6 








290.2 ±0.2 


2.5±0.1 


10 ±53 



Table 5: Expected number of lighter staus to be stopped in the hadronic calorimeter (iV-'^^') and 
the yoke ( N-"^") of the LDC detector at the ILC with ^/s = 600 GeV and an integrated luminosity 
of £ = 250 fb~' for the CMSSM points A D specified in Tabic 3, and expected 68% CL accuracies 
for the determination of m?, Tr and ms for tos — 10 GeV, fa = 10^^ GeV, y = 1 and jegl = 1/3. 
For each CMSSM point, we also list the value of ctsusy, as obtained with Pythia 6.4 [141]. 

the specific shape of the \pT\/mr distribution, which is governed by the mass spectrum and 
the considered ^/s. In particular, for scenario D and ^/s = 600 GeV, almost all the events 
are clustered around the \p~\/mr = 0.2-0.3 region, as can be seen in the bottom right plot 
of Fig. 16. The momenta would thus be so low that the staus do not have enough energy 
to reach the LDC hadronic calorimeter/yoke. Nevertheless, they are energetic enough to 
travel out of the beam pipe and become trapped in the inner detectors. We thus treat all 
the staus produced, 1814, as stopped to derive the expected and accuracies listed 
for scenario D. This is different for the scenarios A-C where we use iVi^^i and Nl°'^^ to 

T T 

infer those accuracies. Here we should note that a future linear collider will have a tunable 
centre of mass energy ^/s. This will allow us to maximize the number of staus stopped in 
those detector parts most suitable to study the decays. However, such optimisations for 
scenarios A-D are beyond the scope of this work and left for future studies. 

Independent of the stopping of staus, the kinematics from the direct production of 
stau pairs will allow for a determination the stau mass with an accuracy listed in Table 5. 
As shown, we expect that the stau mass can be measured very precisely at the ILC, with 
an error below 1% for all four CMSSM points. This can perhaps be further improved on if 
combined with a mass extraction using TOF data. The measurements are not expected 
to be that accurate, but greatly improve with increasing numbers of trapped staus as can be 
seen in Table 5. Accordingly, these measurements would benefit from the ^/s optimisations 
mentioned above. However, the statistical accuracy of the determination expected to 
be achievable at the ILC is already very promising at this stage. 

For scenarios A-C. the explored measurement of the stau lifetime utilises the staus 
trapped both in the hadronic calorimeter and the yoke, while only those trapped in the 
hadronic calorimeter are considered when probing the axino mass. This restriction results 
from the fact that the axino mass exploration relies on the kinematics of the 2-body decay 

TR TO,, 

rria = \J m~ + — 2mrET, (5.1) 

and thereby on the determination of the energy of the emitted tau lepton [37]. With 
Et being the maximum energy of the observable r decay products, the latter requires to 
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study the endpoints of the energy spectrum of the emitted 'r jets' and thus a hadronic 
calorimeter [43]. As shown in Table 5, only rough upper limits on nia are expected for 
rria = 10 GeV. This is due to the difHculties to determine Er with sufficiently high precision. 
While the situation will improve for rria/m^ > 0.1, it will be even worse for smaller ma/m^ 
values because of the demand for an even more precise E^- determination.^^ 

Measurements of the values of t:^, m^f, and mjj can be used to approximately 
constrain the PQ scale fa [37]. Here the mjj dependence becomes only relevant in the 
iT^a/^T > 0.1 region, where one should be able to infer with a better accuracy than 
reported in Table 5. Using the LL expression (3.12), we find that only the following 
combination can be probed by measuring t^, m^, and at colliders 

where one sees explicitly that the dependence becomes negligible for m^/m^^ < 0.1. Note 
that it will not be possible to constrain /„ alone from collider data. For the considered 
SUSY axion models, fa always appears as the ratio Cq/ fa in the prefactor of the decay 
widths and as the product yfa in the logarithm; see (3.2) and (3.4). However, ongoing and 
future axion searches might yield valuable additional insights into the axion model possibly 
realised in nature as detailed, e.g., in Sec. 6 of Ref. [35]. In fact, an axion discovery via the 
axion-photon interaction (2.6) may point to eq and fa- To some extent, ambitious precise 
future measurements of t^, m^, (and mj) at colliders could then probe y via (5.2) and 
thereby - as advertised in the last paragraph of Sec. 3.1 - the hierarchy between /„ and 
the heavy (s)quark masses Mq^q^^. 

5.2 Distinguishing between Axino LSP and Gravitino LSP 

Should a long-lived charged massive particle be detected at a collider, it is not guaranteed 
that its decay products include an axino. Despite being able to accurately determine the 
NLSP's lifetime and other properties, there can remain ambiguities over the identities of its 
non-SM decay products, i.e., the LSP. In particular, as mentioned at the end of Sec. 2, the 
gravitino G is also a well-motivated LSP candidate in SUSY scenarios that include gravity. 
Similarly to the axino, it is very weakly interacting so that a long-lived stau NLSP can also 
occur naturally in gravitino LSP scenarios. Here the stau NLSP lifetime is governed by a 
2-body decay which proceeds via a supergravity coupling [51] 

^GLSP ^ l/r(fR ^ rG) = -^(i--G] , (5.3) 

where the gravitino mass does not only enter via the kinematics but also governs the 
interaction strengths via the spin-1/2 goldstino components. 

Comparing (5.3) with r-^^^ ~ l/r(TR — ^ ra) given by (3.1) in the axino LSP case, 
agreement between the stau NLSP lifetimes in both scenarios is found on the lines shown 
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Figure 17: Agreement between the stau NLSP lifetime in the axino LSP case and the one in the 
gravitino LSP case, t~^^^ = t~^^^, is indicated by the sohd (dashed) hnes for nir = 300 GeV and 
rrig — l.lrrir {'m-^ — 2.0 my). The considered SUSY axion model parameters are jegl — 1/3 (lower 
lines) and {eql — 1 (upper lines) with m~/m~ 1, y — 1, and fa above 6 x 10^ GeV in line with 
the lower limit from astrophysical and cosmological axion studies (1-1). 

in Fig. 17. For concreteness, we consider rri:^ = 300 GeV and = 1.1 (solid lines) 
or 2.0 (dashed lines). Those masses may be known quite precisely from collider studies 
when Tt will be measured. Indeed, the LSP mass nia or will be much harder to probe 
(cf. Table 5) and also fa might still not be more strongly constrained than summarised 
by (1.1). To account for the potential additional eg dependence, the agreement is 
indicated for the two representative cases jegl = 1/3 (lower lines) and jegl = 1 (upper 
lines), whereas y = 1 is kept fixed. In the limit ml ~ <^ m~, we find that 



in (5.3) gives the same stau NLSP lifetime as the LL result (3.12) in the axino LSP case. 

The smaller or fa, the shorter is the stau NLSP lifetime. Accordingly, one finds 
a lower limit imposed by (1.1) in the axino LSP case. In Fig. 17 {m^ = 300 GeV, 
m~/m~ <^ 1), this lower limit ranges between 1.2 x 10"^ s (for = 330 GeV and 
cq = 1/3) and 5.1 x 10^'' s (for rrij^ = 600 GeV and eg = 1). Much smaller r:^ are possible 
with a light gravitino LSP, < 40 keV, and could as such provide hints against the 
axino LSP. In Ref. [37] it had been pointed out that < ms and > days will point to 

i^In a more general setting than the CMSSM, it may also be the e^R that is the NLSP. Then it will be 
easier to achieve the precision required to probe smaller ma values since this relies now on a measurement 
of the energy Ea of the electron that is emitted (instead of the unstable tau) in the decay br — > ea. 




(5.4) 
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a gravitino LSP scenario. While we agree on the short hfetime being a potential indicator 
for the gravitino LSP, the mentioned long lifetime cannot exclude the axino LSP since fa 
may well be larger than 5 x 10^^ GeV considered in [37]. Instead, we find that the upper 
limit < (4-6) x 10^ s = 67-100 min imposed by CBBN for typical yields (4.7) applies 
for both LSP candidates. The additional constraints from hadronic energy injection might 
tighten this upper t:^ limit but differently for the gravitino LSP (see [47,143]) and the axino 
LSP (see [39] and Sec. 4.3 above). 

Another way in which one can try to distinguish between an axino LSP and a gravitino 
LSP is by studying the 3-body decay of the stau tr ^ r + 7 + [37, 44]. In Ref. [37] 
the authors focused on the branching ratio of the integrated 3-body decay rate and on the 
double differential distribution of the 3-body decays. For the axino LSP case, our result 
of the latter quantity is given in (3.4) above and also used to calculate that branching 
ratio (3.7) with results illustrated in Fig. 3. Also Ref. [44] considered the branching ratios 
of the 3-body decays in the following form and with the given choice of cuts, 

r(?R ^ r7X;xf = '-^,xr' = 0.134) 
T{?^^rX) ' ^ = '^'^- 

Here we present an updated comparison of both the double differential distributions 

1 dT{?K^T-fX) 



r(rR ^ r 7 X; x™* = 0.2, xf"^ = 0.134) dxj dcos 6 



X = Z,G, (5.6) 



and of (5.5) using the results of our calculations of the 2- and 3-body decays of the stau 
NLSP into the axino LSP presented in Sec. 3. For the 3-body decay of the stau NLSP into 
the gravitino LSP, we use the result obtained in [37], which is applicable for finite m^M 

Figure 18 shows (5.6) as a function of x^ and cos 9 for ni:^ = 100 GeV and the given 
choice of cuts: x'xj^^ = 0.2 and x^^^ = 0.134. The darker shadings correspond to higher 
number of events, with the contours representing the values 0.2, 0.4, 0.6, 0.8, 1.0, 1.5 and 
2.0. The first three panels show the axino LSP case for = 10^^ GeV, m|/m~ <C 1, 
[cqI = 1/3, y = I and for different values of m^/m^ = 1.01 (top left), 1.1 (top right) and 
2.0 (bottom left). The lower right panel shows the gravitino LSP case for m^/m^ = 2.0 
and niQ = 16 MeV. This value of uiq is chosen such that one finds an agreement of the 
stau lifetime in the two scenarios considered in the bottom panels; cf. (5.4). We however 
note that changes in m^, while significantly altering the stau lifetime, do not greatly affect 
the shown differential distribution as long as m~/m| <C 1. Likewise, changes in the values 
of fa and cause only slight changes of the distributions. 

To varying degrees, all the panels show a preference for producing soft and coUinear 
photons. The top right panel is the scenario most closely resembling the plot shown for the 
axino LSP case in Fig. 5 of Ref. [37]. Here we find an excess of events for highly energetic 
photons produced back-to-back to the emitted taus in comparison to what is obtained in 
the corresponding gravitino LSP case (cf. Fig. 5 of [37]). This is in agreement with the 
result found in [37] , and we now note how this excess region changes for different values of 



Ref. [144] the decay tr — >■ T7G was also considered but in the limit rrig — >■ 00. 
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Figure 18: The normalised double differential distribution (5.6) of the 3-body decays tr — T7X 
with X = a (first three panels) and X = G (bottom right panel) for m^r = 100 GcV. The 
contours shown in the plots represent the values 0.2, 0.4, 0.6, 0.8, 1.0, 1.5 and 2.0, with darker 
shadings corresponding to higher number of events. The cut parameters are set to a;™* = 0.2 
and xg"* = 0.134. For the axino LSP, fa = 10" GeV, m|/m| < 1, jegl = 1/3 and y = 1 in ah 
three panels while mg/mr varies: 1.01 m^; (top left), 1.1 (top right) and 2.0 (bottom left). For the 
gravitino LSP, rrig/m^ = 2.0 and ~ 16 MeV, for which a agreement is encountered in the 
two scenarios considered in the bottom panels; cf. (5.4). 

the bino-stau-mass ratio m^/m^ . For smaller m^/ni:^, the excess is more pronounced with 
events covering a phase space that extends towards larger values of cos 9 as can be seen in 
the top left panel of Fig. 18. This is due to the one-loop decay mode, shown in Fig. 2 (a), 
being more dominant over the two-loop decay mode in that region as the bino becomes more 
and more on-shell. In turn, towards larger values of mj^/m^, the excess region shrinks such 
that for rrij^/m^ = 2.0 (bottom left panel) the excess region is basically missing altogether. 
In the gravitino LSP case, (5.6) is less sensitive to m^/rrir with more events expected in 
the very hard photon region for larger m^/rrir- It must be noted however that this increase 
is very marginal and would require high statistics to be noticeable. Thus, the larger the 
rrij^/m^f ratio, the more difficult it becomes to differentiate between the axino LSP and 
the gravitino LSP. This is shown explicitly in the bottom panels of Fig. 18. Here we have 
chosen parameters for the gravitino LSP case to most closely resemble phenomenologically 
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the axino LSP case to its left to test a particularly difficult situation: nij^/m^f = 2.0 and, 
as already mentioned, such that agrees in the lower panels. If we look closely at 
those two distributions, we find small differences, with the axino LSP scenario having a 
larger phase space area with small values of (5.6). However, these differences are small and 
we would need large statistics and precise measurements to rely on them for a distinction 
between the two scenarios. 

To give an idea of the numbers that can be expected, we assume 5000 trapped staus 
whose decays can be observed in a detector component (such as a hadronic calorimeter) 
that allows for the analysis of the 3-body decays. This may seem somewhat optimistic 
in light of our N-'^''^'^ estimates (obtained without exploring ILC-parameter optimisations) 
given in Table 5. Nevertheless, even iVi^^' = 10^ could be feasible at the ILC already 
without additional stopping-detectors [43] . Also at the LHC one might be able to stop the 
assumed number staus but only when invoking additional stopping-detectors [41,42,44]. 
For example, with massive stoppers set up next to the Compact Muon Solenoid (CMS) 
detector at the LHC, it seems possible to stop staus with ^7 up to 0.6 and to achieve 

^stopped ^^j^gg ^^^^ 

are more than an order of magnitude above the ones reported in 
Table 4 above [44]. Now, for 5000 trapped staus, we expect to see 149±12(stat), 58±8, 
26ib5 and 28ib5 3-body events with 0.2 < x-y < 1 and — 1 < cos 9 < 0.866 when considering 
the scenarios in Fig. 18 in the sequence of left to right and top to bottom, respectively. 
Looking at the region with 0.8 < Xj < 1 and —1 < cos 9 < —0.5, in which an excess 
of events may be found, we would expect 53ib7(stat) (36% of the above number of 3- 
body events), 19±4 (33%), 0.3±0.6 (1.4%), 0.4±0.6 (1.4%) events, respectively. Clearly, 
this looks very pessimistic for distinguishing between the axino LSP and gravitino LSP 
scenario with m^/niTf: = 2.0. If we enlarge the phase space region to — 1 < cos^ < while 
keeping 0.8 < < 1, we would obtain 0.7±0.8 (3.5%) and 1.5ibl.2 (5.2%), respectively, 
for these two cases. Admittedly the improvement is not large, but in lieu of other signals, 
this discrepancy could be indicative. 

In the axino LSP case, the markcdncss of the potential excess of highly energetic 
photons produced back-to-back to the emitted taus is not only sensitive to m^/m^ but 
also to the size of the logarithm log[y^/^/(2m~)]. As can be seen in (3.4), the smaller this 
logarithm, the larger is the relative importance of the one-loop 3-body decay mode shown 
in Fig. 2 (a) and thereby the number of excess events. This can affect the success of the 
distinction between the axino LSP and the gravitino LSP in a way that we now illustrate 
by establishing a connection to the results of Ref. [44]. 

In Ref. [44] the branching ratios (5.5) and distributions in 9 intervals were explored 
based on the 2- and 3-body decay widths calculated in [37]. Accordingly, the results for 
the axino case did rely on a particular choice of ^ log (/a/m), where the factor ^ and the 
mass scale m were introduced to parametrise the uncertainty associated with the cutoff 
procedure applied in [37] instead of the two-loop calculations presented by us in this work; 
see also Sec. 3.2. In particular, in the comparisons of (5.5) in [37,44], log(/a/m) = 20.7 was 
fixed and different choices of including C~yy [37], 0.5, 0.75, 1 and 2 [44] were considered. 
Relying on our 3-body result (3.4) and on the 2-body rate in the LL approximation (3.12), 
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we can now relate those choices to the parameters of the underlying model: 

The different ^ scenarios considered in [37, 44] can thus be related to combinations of y and 
fa for a given m^r. For example, with m^r 2± 100 GeV, ^ ~ 0.5 (1) is found for y = 0.01 
and /a = 6 X 10^ GeV (y = 1 and fa = lO^-'-GeV), while ^ ~ 2 would require y S> 1 
and/or fa 3> 10^^ GeV. With such relations in mind, the results for ^ = 0.5-1 shown in 
Figs. 16-20 of Ref. [44] apply also to our considerations. There it has been estimated, e.g., 
that Ar£*°PP^^ ~ 300 (3000) would be sufficient for ^ = 0.5 (1) to establish a 3a difference 
between the branching ratios (5.5) when 2± 100 GeV and rn^ = 1.1 m^^, while the 
required ]\[p^°PP^'^ increases to about 700 (10^) for = 1.2 m:p. For further results and 
details, we refer to Ref. [44] and to relation (5.7) for the mapping to parameters of the SUSY 
axion models considered in this work. Note however that relation (5.7) is only relevant for 
making contact with the phenomenological investigations of [37,44]. When the logarithm is 
not fixed to a number such as log(/a/m) = 20.7, agreement of our 3-body result (3.4) and 
of the 2-body rate in the LL approximation (3.12) with [37] is found for ^ = 1, m = V^rrir 
and CaYY = Seg, as already discussed in Sec. 3. Still the size of log[y^/a/(2m|)] will be 
relevant for the viability of the axino-gravitino LSP distinction as illustrated by the results 
of Ref. [44] as mentioned above. 

6. Scenarios with mg, < 

In the above sections, we mainly considered the axino LSP scenario with a stau NLSP, 
where < < ruQ. In Sec. 5.2 we also considered the hierarchy < < rria (or 
equivalently, a SUSY scenario with no PQ extension) in the context of differentiating the 
experimental signatures of an axino LSP and a gravitino LSP. Now limiting ourselves to the 
hierarchy nia < rrir, yet another two orderings of the masses are possible, < < 
and uIq < < m^, as already mentioned in Sec. 2. We will explore the consequences of 
both of these in this section. In particular, we consider effects on decays and the lifetime 
of the lighter stau, on cosmological constraints and on collider phenomenology. 

Scenarios in which both the axino and the gravitino are lighter than the other spar- 
ticles (cf. [26, 27]) have been considered first in studies aiming at the understanding of 
structure formation and of the matter power spectrum [145-147]. Moreover, it had been 
realised that such a spectrum may provide a solution to the gravitino problem of thermal 
leptogenesis [77]; see also [63,72,148,149]. Such a setting may also allow for a late non- 
standard increase of ANi, after BBN but before the beginning of structure formation [150]. 
Brief comments on ma,mQ < m^f scenarios were already given in Refs. [35,37]. We now 
elaborate on those and provide several new insights, where our focus is on the behaviour 
of the lighter stau and on associated constraints and phenomenology, see also Ref. [63].^^ 

^^As this work was being finalised, a study of Tr constraints in scenarios with the gravitino LSP and the 
axino NLSP appeared [151], which is complementary to our considerations. 
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Figure 19: The lifetime of the hghter stau (dotted hnes) when it can decay into both the axino 
and the gravitino (ms,TOg < rrir) as a function of mg and /a, for rrir = 300 GeV, rrig = Limy, 
I eg I — 1/3, y = 1 and rn|/m| ^ 1. The sohd grey hnes indicate the relative importance of 
r(TR — 7^ ra) as a fraction x of the total 2-body decay width Tl^f^. Above (below) the x = 0.5 
contour, the stau decays dominantly into the gravitino (axino). 



6.1 Lifetime and Branching Ratios 

Let us first take a look at how the hfetime of the stau is affected. In addition to the 2-body 
decay tr, — >■ ro, the 2-body decay tr — )• tG is also kinematically viable. The lifetime of 
the stau is then governed by these 2-body decays: 



" rlt Tl^'' r(?R^rS) + r(?R^rG) 

with the partial decay widths (3.1) and (5.3). For given m^f and m^^, the relative importance 
of the different 2-body decays depends on fa and m^, which govern r(rR — >• ra) and 
r(?R — >• tG), respectively. This is illustrated in Fig. 19 for m:^ = 300 GeV, = 1.1 m^;, 
legl = 1/3, y = 1 and m~/m~ <C 1, where the latter limit minimises the dependency of 
r(rR — >• ra) on m^. The dotted lines indicate contours of the stau lifetime and the 
solid grey lines the fractional contribution of F-^ = F(rR — )• ra) to the total 2-body decay 
width of the stau: x = F-^/F^q^. ~ BR(rR — t- ra). Note that the upward bending of the x 
contours towards larger results from the factor (1 — m~/rn-~)^ in F(rR — t- tG), cf. (5.3), 
which also accounts for the phase space suppression of the -tr — t- tG decay for — t- m^. 
Above (below) this contour, the stau decays dominantly into the gravitino (axino), which 
explains the clear turning points in the shown lifetime contours. 

For ma,m^ < nir, both 3-body decays -tr — >• Tja/G contribute simultaneously to 
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Figure 20: Branching ratios (dotted lines) of (a) the 3-body decays tr — > rja/G and (b) the 
4-body decays tr — Tqqa/G as a function of mg and fa, for the same parameters as in Fig. 19. 
The cut parameters in (a) are set to a;™* = 20 GeV/m^: ~ 0.07 and Xg^^ — 0.134 and in (b) to 
TO^g* = 2 GeV. The solid grey lines indicate the fractional contributions of the respective stau 
decays into the axino (a) x = T'-' /Tl^f^ and (b) x 



— T^4b 



rf/r: 



T,4b 



tr ^ r + 7 + ^. Accordingly, the corresponding branching ratio with cuts is given by 



BR(3-body;x-*,xr) ^ ^ = ^ " ' %,2b " ' (6.2) 



with (3.6) and the corresponding G result of Ref. [37]. In Fig. 20 (a) this branching ratio 
is illustrated by the dotted lines for the same parameters as in Fig. 19 and for x™* = 
20 GeV /m^f and x™* = 0.134, as in Sec. 5. The shown contours vary only between 0.8% 
and 1.5% and confirm that the total stau decay width F^^^ is governed by the 2-body decays. 
The solid grey lines show the fractional contribution of F-*^ = F(ra — )• T7a; xJ^*^*, x^*^*) to 
the total 3-body decay width of the stau: x = /T\^^ . One finds that F-*^ dominates in 
a (mg, fa) range that is slightly larger than the one in which F^*^ dominates in the 2-body 
case (see Fig. 19). 

Figure 20 (b) shows the corresponding 4-body branching ratio (dotted lines) 

BR(4-body, rugg ) - = , (6.3) 

tot ^ tot 

where _ 

T{Tn^Ta/Gqq;m'gf)= / dm^g- ^ ^ '—^ (6.4) 

Jm'gf drugg 

with (3.11) and the corresponding G result of Ref. [90]. Here we set m.^^* = 2 GeV, as 
in Sec. 4.3, and the other parameters again as in Figs. 19 and 20 (a). Now the solid grey 
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lines show the fractional contribution of = T{t^ — t- raqq^rrf^) to the total 4-body 
decay width of the stau: x = T'-' /Tl^^ . One can sec that T'-' dominates in a {m^, fa) 



region that is even a bit larger than the one in which T'~^ dominates in the 3-body case. 
Moreover, the shown BR(4-body; m^^* = 2 GeV) contours range between 0.04% and 0.16%. 
This demonstrates that the 3-body decays considered in panel (a) are about an order of 
magnitude more abundant than the 4-body decays considered in panel (b). 

6.2 Cosmological Constraints 

Let us now turn to the cosmological constraints from primordial nucleosynthesis for char- 
acteristic freeze-out yields of the lighter stau Yf described by (4.7). Indeed, the freeze-out 
yield of the lighter stau prior to decay is not affected by the presence of a light gravitino. 
Also decays such as G —t- aa or, alternatively, a — t- Ga, even if taking place after freeze OTit 
and prior to decay of the stau, will not dilute as the emitted extremely weakly interact- 
ing particles will not thcrmalise and not affect the entropy density. However, we have just 
seen that can be very different for mj, < m^, and this will affect the CBBN con- 
straints and also the BEN constraints associated with late hadronic and electromagnetic 
energy injection. Moreover, there are now the additional decays tr tG and tr — > rGqq 
contributing to Cem/had- Accordingly, we find 



tot _ 1 
had 

^ tot 



dmqq mqq 

mqq dn^qq 



L 



m^—mpi—mr 



+ / 



dr(?R ^ rGqq) 



^^"^^^ dmqq 



(6.5) 



with r^Qj given in (6.1), and that the conservative estimate (4.30) changes to 

eem = 0.3 BR(?R to) + 0.3 E^ BR(?r tG) , (6.6) 
where the energy of the tau emitted in tr tcl and tr tG is given respectively by 

E": = ^ ^ " ^ and E^ = ^ „ ^ . 6.7 

The resulting constraints are illustrated in the ruQ-fa plane in Fig. 21 for (a) = 
300 GeV and (b) m^f = 1 TeV, with the other parameters set to m~/m!~ <C 1. legl = 1/3, 
y = 1, = 1.1 m^; and the yield = 0.7 x 10~^^ {m:^/l TeV). For orientation, we show 
contours of = 10^, 10"^, and 10^ s (dotted lines) and of x = r|^/r[^f [~ BR(?r ra)] = 
0.01, 0.5 and 0.99 (solid grey lines). In the region with x > 0.99, the (C)BBN constraints 
are insensitive to and basically as in the < < case considered in Sec. 4.3; 
see also Fig. 5 (a) of Ref. [39]. In turn, in the region with x < 0.01, the shown constraints 
are insensitive to and basically as in the G LSP cases with < < or no PQ 
extension; cf. Fig. 6 of Ref. [47]. 

Let us now describe the (C)BBN constraints shown in Fig. 21. The regions to the 
right of and above the long-dash-dotted (red) and the solid lines are disfavoured by the 
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Figure 21: Cosmological constraints on and m^r^ for (a) my = 300 GeV and (b) = 1 TeV and 
rn|/m| <C 1, IbqI = 1/3, ?/ = 1, mj^ = l.lm^ and F?; = 0.7 x 10~^^ TeV). The regions above 

the long-dash-dotted (red) and the solid lines (^Li, ^Be) are disfavoured by the CBBN constraints 
associated with (4.34) and (4.35). No hadronic BBN constraints appear in panel (a). In panel (b) the 
ones associated with (4.26) disfavour the region between the two short-dash-dotted (blue, ^^^^) 
lines. Electromagnetic BBN constraints associated with "^He/D disfavour the regions above the 
double-dash-dotted (green) line. In panel (a) the electromagnetic BBN constraints associated with 
D disfavour the regions between the two leftmost dashed (blue, D^^) lines and to the right of 
the corresponding 3rd line from the left-hand side or the region above the rightmost of those 
lines (Dg™''). In panel (b) the corresponding disfavoured regions are the ones to the right of the 
leftmost dashed (blue, D|^) lines and the ones to the right of the corresponding rightmost (broken, 
DJ;™'') line. Contours of r^; = 10^, 10^, and 10® s are shown by the dotted lines, and contours of 



— -r2b , 



^f/'^lof'l- BR(?R, rS)] = 0.01, 0.5 and 0.99 by the solid grey lines, as in Fig. 19. 



CBBN constraints associated with (4.34) and (4.35). Again, these lines are located close 
to the r;f « 5 X 10^ s contour (not shown). The hadronic BBN constraints associated with 
(4.26) show up in panel (b), where they disfavour the region enclosed in the L-shaped 
strip bounded by the two short-dash-dotted (blue) lines. These severe constraints do not 
appear in panel (a), and also the conservative constraints associated with (4.27) are absent 
in both panels, which is consistent with Fig. 5 (a) of Ref. [39] and Fig. 6 of Ref. [47]. The 
electromagnetic BBN constraints lie well within the region already disfavoured by CBBN in 
our panels (a) and (b). Here the '^He/D constraint disfavours the regions above the double- 
dash-dotted (green) line. Moreover, in panel (a) the D constraint disfavours the regions 
between the two leftmost dashed (blue, D|^) lines and to the right of the corresponding 
3rd line from the left-hand side or the region above the rightmost of those lines (D^^*^), 
whereas in panel (b) the regions to the right of the leftmost dashed (blue, D|^) lines and 
the ones to the right of the corresponding rightmost (broken, Dg^**) line are disfavoured. 

We find that the upper limits on mg, which have been obtained for a given in the 
gravitino LSP case and which imply upper limits on the SUSY breaking scale [47, 143], 
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hold only in the parameter region with x <^ 0.5, i.e., for r(rR — t- ra) <^ r(rR — t- tG). 
For X 3> 0.5, the gravitino mass range 0.1 < < becomes cosmologically viable 
with a standard thermal history, R-parity conservation and a standard freeze-out yield i^. 
Despite a then relatively small branching ratio of the tr tG decay, this might still be 
relevant for collider phenomenology since a kinematic uIq determination appears to be 
feasible only in that mass range [43,44,144] which is disfavoured otherwise [47,143]. On 
the other hand, our upper limits on fa obtained for a given m^f: and [35,39] hold only 
in the parameter region with x 3> 0.5, i.e., for r(TR ra) 3> r(TR — tG). If x <C 0.5, /„ 
values even as large as Mqut — 2 x IQ-*^^ GeV will not be cosmologically disfavoured. 

So we see that the m^, < scenario opens up regions of parameter space that 
are disfavoured if only the axino [35, 39] or only the gravitino [47, 90, 143] is lighter than 
the stau. Accordingly, without additional phcnomcnological insights into or fa from 
collider experiments or axion searches, it will not be possible to infer rigid upper limits 
separately on /„ or niQ once m;^ and are known, even under the assumption of a 
standard thermal history. In fact, for known and and with inferred for a 
standard thermal history with Tr > T[, we think that a figure such as the ones shown in 
Fig. 21 will be the ideal way to display the possible (C)BBN constraints on and/or 
for a given axion model. 

Before proceeding wc also note that some of the constraints disfavour regions for 
"^ai ^ ^^^^ were previously free of constraints. For example, the electromag- 
netic and hadronic BBN constraints in the G LSP cases with mg < < or no PQ 
extension disappear for a finely tuned iTiQ-m^ degeneracy leading to eg^i/had ~^ 0- 
T^ai iT^Q < i^T aiid no additional finely tuned ma-rrir degeneracy, this is no longer the case 
since e^jn/had '^iH be governed by (4.30) and (4.28) which can remain sizeable. Accordingly, 
the D and "^He/D constraints shown in Fig. 21 do not disappear for — > m^. 

6.3 Collider Phenomenology 

Not only the cosmological constraints but also collider phenomenology is affected by the 
possible mass hierarchy m^, < considered in this section. The methods to probe 

and will remain the same. Also our estimates of jv£*°PP'^'l Jsf]}'^^^ and N~°^'^ in Sec. 5.1 
will not change for the assumed T^f values. However, the experimental distinction between 
the axino LSP and the gravitino LSP faces additional complications. While brief comments 
on this issue had already been given in Sees. 1 and 4.3 of [37], we now address this issue 
more extensively and in more detail, with a number of additional insights. 

Even with the axino being the LSP, a short stau lifetime of < 10~^ s is now possible 
for rria < < 40 keV due to r(ra — t- tG) ^ r(rR — t- ra). Then such t-^ will no longer 
provide a hint against the axino LSP (cf. Sec. 5.2) but point to a light gravitino, with 
the possibility of an even lighter axino. Moreover, we have shown explicitly in Fig. 20 (a) 
that both r(TR — >■ raj) S> r(rR — >■ rG'7) and r(TR — )• raj) <C r(rR — )• tG^) is possible, 
whereas the corresponding regions (i.e., x ^> 0.5 and x <C 0.5) in parameter space are 
mainly governed by and fa- Accordingly, double differential distributions such as 
the ones shown in the top panels of Fig. 18 do not necessarily point to the axino LSP 
if m5,m?; < m^r. Instead, they then would point to a parameter region such as the one 



-52- 



with X ^ 0.5 in Fig. 20 (a), on which cosmological constraints can be imposed as shown 
in Fig. 21. The 3-body decays can thereby stih provide very valuable information. Since 
the region with x S> 0.5 in Fig. 20 (a) basicahy coincides with the x S> 0.5 region in 
Fig. 19, an experimental signature for r(TR — > ra'y) S> r(rR tGj) implies r(ra 
TO,) S> r(ra tG) and thereby that the stau decays into the axino govern . This would 
mean that a measurement would probe the PQ scale /a, as outlined in (5.2) for the 
considered SUSY axion models. In turn, the opposite case with x <C 0.5 would mean that 
a measurement would probe the gravitino mass rriQ (and thereby the SUSY breaking 
scale), provided Mpi = 2.4 x 10^^ GeV as inferred from Newton's constant. 

Note that a tt^ measurement alone will point to a tt^ contour in the m^-fa plane; cf. 
dotted lines in Figs. 19 or 21. Such a r^r contour can then be confronted with cosmological 
constraints. However, as can be seen in Fig. 21, the most restrictive and thereby relevant 
constraints basically follow the shape of the contours. Thus, the cosmological constraints 
cannot tell us whether r:^ is governed by stau decays into the axino or the gravitino or 
by both. It will then still be unclear whether a measurement probes /„ or or a 
combination of both. 

Axion searches and resulting limits on /„ may help in particular for small values of tj^. 
Already the current /„ limit (1.1) will point to T(tr tG) » r(rR to) if < 10"'''- 
10^^ s, where the given range of the upper limit reflects an uncertainty with respect to 
the axion model and here the parameters eg and y; cf. Sec. 5.2. In the future, the lower fa 
limit may increase and thereby also the above values of the limit. Alternatively, a future 
axion discovery may point to an /„ value and also to an eq value (cf. Sec. 6 of [35]). Those 
values together with measured m^r, and Tfr can then be tested with relation (5.2). For a 
reasonable y value, which appears only in the logarithm and which will be difficult to probe 
independently of the stau decays (cf. end of Sec. 5.1), an agreement would then point to 
r(rR —7- ra) S> r(TR — t- tG) and thereby to the mg-independent part of the contour; cf. 
dotted lines for x ^ 0.5 in Fig. 19 or 21. In turn, a severe mismatch in (5.2) may originate 
from r(TR — )• tG) ^ r(TR, — )■ ra) for which is governed by and independent of /„; 
cf. dotted lines for x ^ 0.5 in Figs. 19 or 21. 

Additional information might be gleaned from the kinematics of the 2-body decays 
tr rajG. If Txia and are sufficiently different, two isolated peaks may be observed in 
the energy spectrum of the tau, as given by (6.7). The heights of the peaks are goverend 
by r(TR — 5- to) and r(TR — t- tG). This information, together with (5.3), could in principle 
allow us to determine which of the two peaks corresponds to the gravitino, and thus whether 
the axino or the gravitino is the LSP. This procedure will not work, however, if one of the 
decay channels is strongly suppressed. 

While this may sound encouraging, one faces the complication of the short tau lifetime 
and that the decay products of the tau always include neutrinos. Thus, the two E^^'^ 
values will not show up in the form of two basically monochromatic peaks but, e.g., as tau 
jet energy distributions, which will be further smeared out by the finite detector energy 
resolution as shown, for example, in Fig. 8 of Ref. [44]. Because of such complications, it 
had been pointed out before that kinematical determinations of [37] and [144] will 
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be feasible only for 0.1 m:^ < "^a/G ^ ^t"' However, maybe it is not the stau tr that is the 
lightest MSSM superparticle but the right-handed selectron cr, which would then decay 
into an electron instead of the tau: br — > ea/G. This would facilitate the above also for 
smaller m~^Q since the corresponding two discrete values of Ee^^ would be deposited fully 
in the calorimeters. In the absence of a degeneracy of the axino and gravitino masses, one 
may then observe the two decay modes very distinctively and thereby identify the LSP and 
the decay(s) that govern r^f. 

7. Conclusions 

Supersymmetric Peccei-Quinn models, where both the strong CP problem and the hierar- 
chy problem are solved simultancoTisly, arc compelling also from the cosmological point of 
view. They include attractive candidates for dark matter. On the one hand, dark matter 
can just reside in the form of axions [6-10]. On the other hand, the axino a can be the 
stable LSP in an R-parity-conserving realisation of SUSY and as such be another natural 
dark matter candidate [22,28-36]. In this work we have studied in detail the cosmological 
effects and potential collider phenomenology of such a model, in which the axino is the 
LSP and a charged slepton/stau the NLSP. 

Wc focus on hadronic KSVZ axion models [11,12,14,15], where the axion supcrmultiplet 
couples to the MSSM fields only via loops of very heavy (s)quarks, whose mass is of the 
order of the breaking scale /„ of the PQ symmetry. In Sec. 3 we have calculated the 2- 
body decay tr ra, which is induced only at the two-loop level. Here the large mass 
hierarchy between the KSVZ fields and the MSSM fields has allowed us to use the heavy 
mass expansion technique. Due to a large logarithm, the rate of this decay mode is large 
and dominates all other modes by a factor of 100 or so. We have also calculated differential 
rates for the 3-body decay tr -> raj, which is relevant for collider prospects, and 4-body 
decays such as tr — >■ raqq, which are relevant for BBN constraints. Explicit expressions 
are given for rates and distributions, and their dependencies on the involved parameters 
have been illustrated. For the natural mass range of the stau NLSP (and the other MSSM 
sparticles), its lifetime can be typically of the order of 0.01 s to thousands of seconds. 

In Sec. 4 we have aimed at a detailed and complete account of the cosmological con- 
straints on axino dark matter scenarios with such a long-lived charged slepton. Covering 
a wide range of the axino mass and of the reheating temperature after inflation Tr, 
various possible axino populations have been considered: thermal relic axinos, thermally 
produced axinos and non-thermally produced axinos from NLSP decays. We have ex- 
plored the constraints on those populations imposed by the present dark matter density, 
by structure formation (and associated limits on the corresponding present free-streaming 
velocity) and by bounds on non-standard contributions to the radiation density, AA^^,, at 
the onset of BBN. At present, we find that the BBN bounds on AiV^ are such that they 
do not impose relevant restrictions on the considered scenarios. This is different for the 
dark matter density which leads to upper limits on Tr, and m^. The detailed form of 
these limits depends significantly on /„ and also on the other parameters as explored and 
discussed in Sec. 4. 
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Other BBN constraints relate to the primordial stau NLSP abundance after freeze out 
and prior to decay. Here we have studied various scenarios with and without stau-slepton 
coannihilation, stau-bino coannihilation and simultaneous stau-bino-slepton coannihila- 
tion, leading to different values of the stau yield 1?^. For /„ > 10^^ GeV, the stau decays so 
late that its decay products can reprocess BBN or that negatively charged staus can form 
bound states leading to catalysed BBN (CBBN) of ^Li and ^Be. We have shown explicitly 
that the resulting constraints from late hadronic energy injection and from CBBN are then 
often the most relevant ones disfavouring significant regions in parameter space. 

Those BBN constraints also apply to axino hot dark matter scenarios with a charged 
slepton NLSP and an axion condensate providing the cold dark matter [33]. For such 
light hot thermal relic axinos, structure formation furthermore imposes the limit < 
37 eV. Assuming instead that (almost) all dark matter resides in the form of non-thermally 
produced axinos, the constraints imposed by structure formation considerations can even 
disfavour mj smaller than 1-100 GeV. 

The precise position of the cosmological constraints depends not only on /„, and 
the other parameters but also on the adopted limits on the free-streaming velocity or the 
abundances of primordial nuclei. Since both types of limits rely on astrophysical studies, 
which are still subject to ongoing research, we have considered conservative and restrictive 
limits when presenting the constraints from structure formation and BBN. Taking even 
the most restrictive limits into account, we find that axino dark matter scenarios with a 
long-lived stau remain a viable possibility that can be compatible with the cosmological 
constraints, in particular, for /„ < 10^^ GeV and m:p within the discovery reach of the LHC 
and a future linear collider such as the ILC. 

In Sec. 5 we have explored the prospects to study scenarios with the axino LSP and 
a stau NLSP in collider experiments. Evaluating the stau decay length in the axino LSP 
case, we find that it is typically much larger than the detector size. While most staus (if 
produced) will thus decay outside of the detectors, a fraction of slow staus can get stopped 
in the main detectors or in some additional dedicated stopping material. Those stopped 
staus could then allow us to study their decays to determine their lifetime and the axino 
mass TTT-a, to probe fa and to distinguish between axino LSP and gravitino LSP. 

To quantify the collider prospects in a representative way, we have studied four bench- 
mark scenarios with CMSSM spectra that cover the aforementioned coannihilation scenar- 
ios. We have calculated the stau production rates and the number of staus expected to get 
stopped in a typical detector at the LHC and the ILC. We find that with moderate amounts 
of luminosity, several hundred staus may be stopped at the LHC, which would allow one to 
determine their lifetime with about 10% uncertainty. At the ILC, these numbers improve 
to several thousand stopped staus and an estimated uncertainty of about 3%, respectively. 
Higher statistics are conceivable with additional dedicated stopping detectors at the LHC 
or by tuning the beam energy at the ILC. 

While the stau lifetime and mass can be measured fairly precisely, the determination 
of the properties of the LSP that it decays into is more challenging. We have addressed 
the question whether it is possible to experimentally distinguish an axino LSP from an 
alternative scenario with a gravitino LSP. A short decay time T^f < 10~ ''-10"^ s may point 
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toward the gravitino LSP case, but in general a t:^ measurement alone cannot clearly 
discriminate between an axino LSP and a gravitino LSP. On the other hand, the 3-body 
decays tr T'ya/G can help to distinguish between these cases for large parts of the 
parameter space. However, the distinction will also become challenging in this channel for 
a heavy bino or large values of the ratio jalfn^- 

Finally, we have studied scenarios in which both the axino and the gravitino are lighter 
than the charged slcpton/stau that is assumed to be the lightest MSSM sparticlc. Such a 
mass hierarchy is well motivated and may occur, e.g., in models with gauge- mediated SUSY 
breaking [26,27]. We have shown explicitly the way in which the simultaneous existence 
of a light axino and a light gravitino affects the 2-body, 3-body and 4-body stau decays 
and thereby the stau lifetime, cosmological constraints and collider prospects. Our findings 
show that can now be much smaller and the cosmological constraints can be much weaker 
than expected without a second extremely weakly interacting sparticle being lighter than 
the MSSM sparticles. In particular, values even as large as Mgut 2 x 10^^ GeV 
- which are disfavoured for < < mg [39] - can now be allowed for mg < 10- 
100 GeV < m^. Alternatively, the gravitino mass range 0.1 < iuq < can become 
cosmologically viable with a standard thermal history and a standard freeze-out yield 
if fa < 10^^ GeV and < m^r. 

For a scenario with m^, mg < m^, the stau lifetime depends on the PQ scale fa as 
well as on the gravitino mass mg (and the Planck mass Mpi). Thus a measurement of 
yields only ambiguous information. We have shown that this ambiguity may be resolved 
by extracting the branching ratios of decays into axinos and into gravitinos. Here we 
have proposed two complementary approaches one relying on a study of the 3-body decays 
tr — >■ T^a/G and the other on a study of the 2-body decay kinematics. In fact, the outlined 
kinematical considerations may clarify also the mass hierarchy between the axino and the 
gravitino. In practice this requires the identification of a double-peak feature in the tau 
energy spectrum, which however will only be feasible if \mg — m^l is relatively large. Such 
an analysis would become much simpler with electrons in the case that the NLSP is a 
selectron instead of a stau. It remains to be seen how far detector technology can be 
driven in order to study the various aspects of late decays of charged sleptons. 

In the present exciting time in which the LHC experiments are taking data, there are 
high hopes for nearby discoveries of new physics and new particles. In fact, maybe a long- 
lived charged massive particle will be seen soon in the detectors. Then our detailed study 
presented in this work may become very relevant. This will be even more so if in addition 
also the axion is discovered in ongoing or future axion searches. It may then very well be 
axions and axinos that make up the dark matter in our universe. 
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A. Heavy Mass Expansion and Reduction to Basic Two-Loop Integrals 



A.l Heavy Mass Expansion 

Heavy mass expansion (HME) is an asymptotic expansion in small masses and momenta 
[152, 153]. It applies in the cases where the following assumptions can be made 

• each mass can be categorised as either large, Mj, or small, mf, 

• all external momenta are small compared to the large masses, ki <C Afj. 

For a Feynman integral T of Feynman diagram, F the expansion is compactly written as 



where 7 represents the subgraphs over which the sum is performed (a subgraph is derived 
by drawing all possible one-particle irreducible graphs containing the heavy propagators). 
The operator T performs a Taylor expansion on the subgraphs with respect to the small 
parameters. J'^(r\7) represents the original graphs where the subgraphs have been shrunk 
to a point. 

When considering processes involving loops, we must pay special care to the loop 
momenta qi. Since these are unconstrained by momentum and energy conservation, they 
can be both of the order of the small parameters and of the large masses. In order to carry 
out the expansion in a consistent manner, we must look at each of the momenta regions, 
individually. 

For the two-loop processes considered in this paper, there are four such regions. Let 
us refer to the diagrams in Fig. 1 and denote the momentum running through the outer 
loop qi, the momentum running through the heavy loop q2 and the heavy mass scale M. 
Then these 4 regions are defined as follows: 

1. \qi\,\q2\ < M 

2. -M, |g2| <M 

3. \qi\ < M, 1^2! ~ M 
4- \qi\, \q2\ ^ M 

After the appropriate Taylor expansions, the first two regions reduce to scaleless integrals, 
which vanish in dimensional regularisation. This leaves only the third and fourth regions, 
which in practical terms translate into the Taylor expansion of only the heavy loop and the 
Taylor expansion of the whole diagram with respect to the small parameters, respectively. 
Figure 22 gives a pictorial depiction of what is going on, where the propagators with the 
large mass parameter M are represented by thick lines and the rest by thin lines. 




(A.1) 



7 
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Figure 22: The leading term in tlie heavy mass expansion of the two-loop 2-body decay of the stau 
NLSP into the axino LSP, where the propagators with the large mass parameter M are represented 
by thick lines and the rest by thin lines. 



A. 2 Tensor Reduction 

Using Passarino-Veltman reduction procedures, the tensor integrals were simplified to 
scalar ones. Using the following notation, 



ToK,n„n3) = J _M|^(,i -M|)-fe +g,)-3 ' (^-2) 

we make use of the following relations in order to reduce our tensor integrals: 

r^^'^nm,^2,n3) = ^9^''''qi.q2To{nun2.n^), (A.4) 

+{D{qi-q2f - qlqlW'g^'''' + g^''''9'''''')\To{ni,n2^n^), (A.5) 
where g^^'^ is the metric tensor and D is the number of dimensions. 

A. 3 Integration by Part Identities 

After applying the heavy mass expansion and tensor reduction formulae, we are left only 
with vacuum two-loop integrals and one- loop integrals (the results of which are well-known) . 
The two-loop integrals take the form To(?t.i, n2, na). These are simplified with the use of 
integration by parts (IBP) identities, reducing the two- loop integrals to master integrals. 
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Using the following identities, 

d^qi d'^q2 TTU , ,n , , o ^\ , , ^ = (A.6) 

d^9i rf^92 A ,., „2 'l.2^..... , = (A-9) 



we reduced all our integrals to master integrals. This was implemented using the Auto- 
mated Integral Reduction program [154]. In fact, in the course of our calculation, we only 
had one master integral, ro(l, 1,0), where the two-loop integral decouples into a product 
of two one-loop integrals. 
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